ABSTRACT Title of Document: A MATHEMATICAL FRAMEWORK FOR OPTIMIZING DISASTER RELIEF LOGISTICS Abbas Mohasel Afshar, 2011 Directed By: Professor Ali Haghani Department of Civil and Environmental Engineering In today?s society that disasters seem to be striking all corners of the globe, the importance of emergency management is undeniable. Much human loss and unnecessary destruction of infrastructure can be avoided with better planning and foresight. When a disaster strikes, various aid organizations often face significant problems of transporting large amounts of many different commodities including food, clothing, medicine, medical supplies, machinery, and personnel from several points of origin to a number of destinations in the disaster areas. The transportation of supplies and relief personnel must be done quickly and efficiently to maximize the survival rate of the affected population. The goal of this research is to develop a comprehensive model that describes the integrated logistics operations in response to natural disasters at the operational level. The proposed mathematical model integrates three main components. First, it controls the flow of several relief commodities from sources through the supply chain until they are delivered to the hands of recipients. Second, it considers a large-scale unconventional vehicle routing problem with mixed pickup and delivery schedules for multiple transportation modes. And third, following FEMA?s complex logistics structure, a special facility location problem is considered that involves four layers of temporary facilities at the federal and state levels. Such integrated model provides the opportunity for a centralized operation plan that can effectively eliminate delays and assign the limited resources in a way that is optimal for the entire system. The proposed model is a large-scale mixed integer program. To solve the model, two sets of heuristic algorithms are proposed. For solving the multi- echelon facility location problem, four heuristic approaches are proposed. Also four heuristic algorithms are proposed to solve the general integer vehicle routing problem. Overall, the proposed heuristics could efficiently find optimal or near optimal solution in minutes of CPU time where solving the same problems with a commercial solver needed hours of computation time. Numerical case studies and extensive sensitivity analysis are conducted to evaluate the properties of the model and solution algorithms. The numerical analysis indicated the capabilities of the model to handle large-scale relief operations with adequate details. Solution algorithms were tested for several random generated cases and showed robustness in solution quality as well as computation time. A MATHEMATICAL FRAMEWORK FOR OPTIMIZING DISASTER RELIEF LOGISTICS By Abbas Mohasel Afshar Dissertation proposal submitted to the Faculty of the Graduate School of the University of Maryland, College Park, in partial fulfillment of the requirements for the degree of Doctor of Philosophy 2011 Advisory Committee: Professor Ali Haghani, Chair Professor Gerald Galloway Professor Bruce Golden Professor Paul Schonfeld Professor Lei Zhang ? Copyright by ABBAS MOHASEL AFSHAR 2011 ii Dedication I dedicate this dissertation to all the wonderful people in my life? Especially to my parents for their unconditional love, To my lovely wife for her support and patience, To my teachers for instilling the power of higher education within me, To my friends for supporting me through this journey, And to my daughter to be who has inspired me to work harder than ever. iii Acknowledgements I would like to sincerely express my gratitude to my advisor professor Ali Haghani for his continued support and assistance throughout my Ph.D. studies in the University of Maryland. Special thanks to the members of the supervisory committee, Professor Gerald Galloway, Professor Bruce Golden, Professor Paul Schonfeld and Professor Lei Zhang for their guidance and support during the completion of this dissertation. I offer my regards and blessings to my student colleagues who helped me in so many aspects during the course of my studies. This dissertation would never have been completed without the encouragement and assistance of my beloved family and friends. iv Table of Contents Dedication ............................................................................................................... ii Acknowledgements ................................................................................................ iii Table of Contents ................................................................................................... iv List of Tables ........................................................................................................ vii List of Figures ...................................................................................................... viii List of Abbreviations ...............................................................................................x Chapter 1: Introduction ............................................................................................1 1.1 Disasters .................................................................................................1 1.1.1 Definitions.....................................................................................1 1.1.2 Numbers and Trends .....................................................................4 1.2 Emergency Management .......................................................................6 1.3 Federal Emergency Management Agency .............................................9 1.4 FEMA?s Logistics Supply Chain .........................................................14 1.5 Motivation and Objective of the Research ...........................................18 1.6 Contributions of the research ...............................................................19 1.7 Organization of the dissertation ...........................................................21 Chapter 2: Literature review ..................................................................................22 2.1 Supply Chain Management ..................................................................22 2.2 Facility Location Problem....................................................................26 2.2.1 P-Median Problem ......................................................................27 2.2.2 Covering Problem .......................................................................29 2.2.3 P-Center Problem ........................................................................32 2.3 Vehicle Routing Problem .....................................................................33 2.3.1 VRP Mathematical Formulation .................................................35 2.3.2 VRP Variants ..............................................................................37 2.4 Commercial Supply Chain vs. Emergency Response Logistics ..........38 2.5 Logistics in Disaster Response ............................................................42 2.5.1 Early Ages ...................................................................................43 2.5.2 Recent Studies .............................................................................44 2.6 Conclusions ..........................................................................................47 Chapter 3: Problem Description and Formulation .................................................49 3.1 Problem Description ............................................................................49 3.2 Time-Space Network ...........................................................................54 3.3. Modeling Approach ............................................................................56 3.4 Assumptions .........................................................................................58 3.5 Mathematical model.............................................................................60 3.5.1 Notations .....................................................................................60 3.5.2 Parameters ...................................................................................61 3.5.3 Decision Variables ......................................................................62 v 3.5.4 Objective Function ......................................................................63 3.5.5 Constraints ..................................................................................63 3.6. Summary .............................................................................................68 Chapter 4: Preliminary Numerical Study ...............................................................69 4.1 Design of the sample problems ............................................................69 4.1.1 NETWORK.................................................................................71 4.1.2 Supply and Demand ....................................................................74 4.1.3 Transportation .............................................................................76 4.1.4 Network links and travel times ...................................................76 4.1.5 Time Scale ..................................................................................77 4.1.6 Geographical Decomposition......................................................79 4.2 Generating formulation for commercial solver ....................................81 4.3 Numerical results and analysis.............................................................82 4.4. Summary of the preliminary numerical experiments..........................93 Chapter 5: Solution Approaches ............................................................................94 5.1 General Solution Approaches for Integer Programs ............................94 5.1.1 Cutting Planes .............................................................................95 5.1.2 Branch and Bound.......................................................................96 5.1.3 Lagrangian Relaxation ................................................................98 5.1.4 Benders Decomposition ..............................................................99 5.1.5 Heuristics ..................................................................................101 5.1.6 Structural Exploitation ..............................................................102 5.2 Solution approaches used in emergency logistics literature ..............103 5.3 Solution Techniques for Proposed Mathematical Model ...................108 5.4 Algorithms for solving the location problem .....................................109 5.4.1 Explicit Enumeration ................................................................110 5.4.2 Branch and bound - Hierarchical decomposition......................111 5.4.3 Highest Capacity Ratio .............................................................114 5.4.4 Static Network-Location ...........................................................116 5.5 Algorithms for solving Vehicle Routing Problem .............................120 5.5.1 T-Counter Heuristic ..................................................................120 5.5.2 Origin-Based T-Counter Heuristic ............................................122 5.5.3 Y-List Heuristic ........................................................................123 5.5.4 Y-List Modal Heuristic .............................................................125 5.5.5 Comparing Performance of the Proposed Algorithms ..............126 5.6 Testing robustness of the proposed VRP heuristics ...........................132 5.7 Summary ............................................................................................134 Chapter 6: Detailed Analysis of the Mathematical Model...................................137 6.1 Introduction of the numerical case study ...........................................137 6.2 Sensitivity Analysis of the Structural Parameters ..............................139 6.2.1 Analysis of Number of Commodities .......................................139 vi 6.2.2 Analysis of the Number of Transportation Modes ...................141 6.2.3 Analysis of Time-Step t ............................................................146 6.3 Sensitivity Analysis of the Main Input Parameters............................150 6.3.1 Sensitivity Analysis for Parameters of the Facility Location ...151 6.3.2 Sensitivity analysis for Parameters of Vehicle Routing ...........161 6.3.3 Sensitivity analysis for Parameters of the Commodity Flow....168 6.4 Summary and Conclusions ................................................................178 Chapter 7: Prepositioning and Equity ..................................................................184 7.1 Prepositioning of supplies and vehicles .............................................184 7.1.1 Mathematical Formulation ........................................................185 7.1.2 Numerical Experiments ............................................................187 7.2 Equity .................................................................................................192 7.2.1 New Equity Constraints ............................................................194 7.2.2 Numerical Experiment ..............................................................198 Chapter 8: Summary and Future Research ..........................................................203 Recommendations for Future Research ...................................................215 Bibliography ........................................................................................................220 vii List of Tables Table 1.1 Emergency support function annexes of the national response framework Table 2.1 Commercial supply chains vs. humanitarian relief chains Table 4.1 List of facilities in distribution network Table 4.2 List of required items for survivors of a disaster Table 4.3 Descriptions of the numerical case studies Table 4.4 Summary of optimization results for preliminary experiment Table 5.1 Branch and bound and hierarchical decomposition Table 5.2 Values of Loci variables for HCR heuristic algorithm Table 5.3 Model statistics and optimization results from CPLEX solver Table 5.4 Numerical results of the proposed VRP heuristic algorithms Table 5.5 Comparative rankings of the proposed heuristic algorithms Table 5.6 Objective function values for ten random cases Table 5.7 Running time comparisons of the algorithms for ten random cases Table 6.1 Problem sizes for different number of commodities Table 6.2 Problem sizes for different number of transportation modes Table 6.4 Problem sizes for different values of time step t Table 6.5 Optimization results for different time steps Table 6.6 Analysis of the number of the mobilization centers Table 6.7 Analysis of the number of the federal operational staging areas Table 6.8 Analysis of the number of the state staging areas Table 6.9 Sensitivity analysis on the total number of temporary facilities Table 6.10 Loading and unloading capacities Table 6.11 Optimization results for various loading and unloading capacities Table 6.12 Sensitivity analysis for the number of vehicles Table 6.13 Sensitivity analysis for the vehicle's capacity Table 6.14 Sensitivity analysis for the network travel speeds Table 6.15 Sensitivity analysis for supply availability Table 6.16 Sensitivity analysis for variations in demand Table 6.17 Sensitivity analysis for the relative urgency factor Table 7.1 Introduction of the prepositioning case studies Table 7.2 Percent of supply prepositioned at each facility type Table 7.3 Number of vehicles prepositioned at each facility type Table 7.4 Effect of equity constraint on demand satisfaction at PODs Table 7.5 Relative unsatisfied demand for various equity tolerance levels viii List of Figures Figure 1.1 Number of the reported natural disasters per year around the world Figure 1.2 Number of the victims of the natural disasters per year Figure 1.3 Economic damage of the natural disasters over time Figure 1.4 Four phases of the emergency management cycle Figure 1.5 FEMA?s supply chain structure Figure 2.1 Supply chain process Figure 3.1 A sample physical network Figure 3.2 A sample time-space network Figure 4.1 General disaster area of the numerical study Figure 4.2 Disaster areas in the two affected states Figure 4.3 Map of the federal level facilities Figure 4.4 Map of the state level facilities in Mississippi Figure 4.5 Map of the state level facilities in Louisiana Figure 4.6 The issue of scale in the disaster area Figure 4.7 Geographical decomposition to implement two time steps Figure 4.8 Connecting two time-space networks with different time steps Figure 4.9 Percent of unsatisfied demand over time for case 1 Figure 4.10 Percent of unsatisfied demand over time for case 2 Figure 4.11 Percent of unsatisfied demand over time for case 3 Figure 4.12 Percent of unsatisfied demand over time for case 4 Figure 4.13 Percent of unsatisfied demand over time for case 5 Figure 4.14 Percent of unsatisfied demand over time for case 7 Figure 4.15 Percent of unsatisfied demand over time for case 8 Figure 4.16 Percent of unsatisfied demand over time for case 9 Figure 4.17 Percent of unsatisfied demand over time for case 10 Figure 5.1 Convergence rate comparison of the proposed algorithms Figure 6.1 Disaster area map and facility locations Figure 6.2 Problem size and solution time versus number of commodities Figure 6.3 Comparison of performance for single and multi-modal cases Figure 6.4 Variations of the problem size and CPU time for different time steps Figure 6.5 Effect of the number of the mobilization centers on the objective function value Figure 6.6 Effect of the number of the federal operational staging areas on objective function value Figure 6.7 Effect of the number of the state staging areas on the objective function Figure 6.8 Objective function and CPU time versus total number of temporary facilities ix Figure 6.9 Objective function values versus variations of the loading and unloading capacities Figure 6.10 Effect of the number of vehicles on the objective function and CPU time Figure 6.11 Sensitivity analysis for the vehicle's capacity Figure 6.12 Sensitivity analysis for the amount of available supplies Figure 6.13 Optimization results for random demand values Figure 6.14 Sensitivity analysis of the relative urgency factor for each commodity Figure 6.15 Effect of relative urgency factor on variations of the unsatisfied demand over time Figure 6.16 Effect of using high priority for POD 29 Figure 6.17 Effect of using high priority for POD 41 Figure 7.1 Satisfaction rate variations among points of distribution x List of Abbreviations CRED Center for Research on the Epidemiology of Disasters CSS Commercial Storage Site DHS Department of Homeland Security DLA Defense Logistics Agency ESF Emergency Support Function FEMA Federal Emergency Management Agency FOSA Federal Operational Staging Area GSA General Services Administration LC Logistics Center MOB Mobilization centers NRF National Response Framework PAHO Pan American Health Organization POD Point of Distribution SSA State Staging Area WHO World Health Organization 1 Chapter 1: Introduction In this chapter a general introduction of disasters and disaster management concepts is presented. Section 1.1 provides some general definitions of disasters and some disaster numbers and trends in recent years. Section 1.2 introduces emergency management. Then in section 1.3 federal emergency management agency (FEMA) is introduced. FEMA?s logistic supply chain is discussed in section 1.4. Motivation and objective of this research is emphasized in section 1.5 followed by the contributions of this research in section 1.6. Finally, the organization of the rest of this dissertation is summarized in section 1.7. 1.1 Disasters In this section, first definitions of disasters are presented followed by some statistics from large-scale disasters in the recent years. 1.1.1 Definitions The term ?disaster? is usually applied to a breakdown in the normal functioning of a community that has a significant adverse impact on people, their works, and their environment, and overwhelms local response capacity. This situation may be the result of a natural event such as a hurricane or earthquake; or it may be the result of human activities (PAHO 2001). Some organizations make a distinction between ?disasters??the result of natural phenomena?and ?complex emergencies? that are the product of armed conflicts or large-scale violence and often lead to massive displacements of people, famine, and outflows of refugees. 2 A disaster, as defined by the World Health Organization (WHO), is any occurrence that causes damage, destruction, ecological disruption, loss of human life, human suffering, deterioration of health and health services on a scale sufficient to warrant an extraordinary response from outside the affected community or area. The American Red Cross defines a disaster as an occurrence or situation that causes human suffering or creates human needs that the victims cannot alleviate without assistance. Earthquakes, hurricanes, tornadoes, volcanic eruptions, wild fires, floods, blizzard, drought, terrorism, chemical spills and nuclear accidents are included among the causes of disasters, and all have significant devastating effects in terms of human injuries and property damage. Alexander (1999) defines natural disaster as some rapid, instantaneous or profound impact of the natural environment upon the socio-economic system. He also recommends Turner?s (1976) definition of natural disaster as ?an event, concentrated in time and space, which threatens a society or subdivision of a society with major unwanted consequences as a result of the collapse of precautions which had previously been culturally accepted as adequate?. Center for Research on the Epidemiology of Disasters (CRED), collaborating center with WHO and United Nations, defines disaster as ?A situation or event, which overwhelms local capacity, necessitating a request to national or international level for external assistance; an unforeseen and often sudden event that causes great damage, destruction and human suffering? (CRED 2007). 3 The official definition of disasters in the United States is presented in the Stafford Act. The Robert T. Stafford Disaster Relief and Emergency Assistance Act is the primary legislation in the United States authorizing the federal government to provide disaster assistance to states, local governments, families, and individuals. The Stafford Act defines a disaster as ?Any natural catastrophe (including hurricane, tornado, storm, high water, wind driven water, tidal wave, tsunami, earthquake, volcanic eruption, landslide, mudslide, snowstorm or drought), or, regardless of cause, any fire, flood or explosion, in any part of the United States, which in the determination of the President causes damage of sufficient severity and magnitude to warrant major disaster assistance under this Act to supplement the efforts and available resources of States, local governments, and disaster relief organizations, in alleviating the damage, loss, hardship, or suffering caused thereby.? As these definitions indicate, a disaster is a ?catastrophe? of such magnitude and severity that the capacities of states and local governments are overwhelmed. So the threshold for determining what constitutes a disaster depends upon the availability of resources and capabilities of responding communities. Consequently, a disaster can be prevented by increasing the capacity of the responding organizations. 4 1.1.2 Numbers and Trends From global perspective, the number of natural disasters is increasing every year. For example in 2005, there was 489 country-level disasters affecting 127 countries around the globe resulting in 104,698 people killed and 160 million people affected. For the same year of 2005, the economic damage estimate varies from 159 billion to 210 billion in US dollars. Because of the population growth and new developments in risk prone regions, the exposure of the human kind to the natural disasters is increasing even more. Figure 1.1 shows the number of reported natural disasters around the globe from 1980 to 2007. A least-square linear regression trend-line is drawn to better illustrate the overall pattern. Trend-line in Figure 1.1 shows that in spite of fluctuations due to cyclic or seasonal patterns, the average number of disasters is growing in the long term. During 1980s number of disasters was around 180 per year on average. In 1990s, the average number of disasters was increased to around 300 per year. And in the 2000-2007 period, it was around 460 disasters per year which indicates a dramatic increase. An increase of this magnitude can be explained partially by the global warming theory, and partially by the attention of the media which has increased the numbers of reported disasters all over the world. As the number of disasters increases every year, more people are affected by these disasters. Figure 1.2 illustrates the number of victims of natural or man- made disasters in the last twenty years. The number of victims includes the people killed, injured, lost their homes or evacuated as a direct result of the disaster. As 5 can be seen in figure 1.2, the number of victims has higher fluctuations over the years. However, the trend-line shows a slow increase in the average number of peoples affected each year over time. The number of victims is generally between 100 million and 400 million per year. The exceptionally high number in 2002 is due to a drought solely affecting 360 million in India and China and a major wind storm and flood affecting 160 million people in China. Figure 1.1 Number of reported natural disasters per year around the world (CRED) Figure 1.2 Number of victims of natural disasters per year (CRED) 0 100 200 300 400 500 600 1980 1983 1986 1989 1992 1995 1998 2001 2004 2007 Di sas ter s per ye ar 0 100 200 300 400 500 600 700 1980 1983 1986 1989 1992 1995 1998 2001 2004 2007 Per son A ffec ted (M illi on) 6 Another important factor is the monetary cost of the natural disasters. Figure 1.3 shows the amounts of global economical damage caused by the natural disasters from 1980 to 2007. The average cost per year is $45 billion from 1980 to 1999. However, for 2000 to 2007 period, the average cost is more than $80 billion per year. The linear trend-line shows an increase in the economical damage of the natural disasters over time. Two major disasters affecting the trend are the Kobe earthquake in 1995 and hurricane Katrina in 2005. Figure 1.3 Economic damage of the natural disasters over time (CRED) 1.2 Emergency Management Emergency management (or disaster management) is the discipline of avoiding risks and dealing with risks (Haddow et al. 2007). No country and no community are immune from the risk of the disasters. However, it is possible to prepare for, respond to and recover from disasters and limit the destructions to a certain degree. Emergency management is a discipline that involves preparing for disasters before they happen, responding to disasters immediately, as well as 0 50 100 150 200 250 1980 1983 1986 1989 1992 1995 1998 2001 2004 2007 Billi on US D oll ars 7 supporting, and rebuilding societies after the natural or man-made disasters have occurred. Emergency management is a continuous process. It is essential to have comprehensive emergency plans and evaluate and improve the plans continuously. The related activities are usually classified as four phases of Preparedness, Response, Recovery, and Mitigation. Figure 1.4 illustrates the order of these phases according to the onset of the disaster. Appropriate actions at all points in the cycle lead to greater preparedness, better warnings, reduced vulnerability or the prevention of disasters during the next iteration of the cycle. Figure 1.4 Four Phases of Emergency Management Cycle Some of the main activities during the four phases of the emergency management cycle are summarized below: Preparedness ? Activities to improve the ability to respond quickly in the immediate aftermath of an incident. ? Development of response procedures, design and installation of warning systems, evacuation planning, exercises to test emergency operations, and training of emergency personnel. 8 Response ? Activities during or immediately following a disaster to meet the urgent needs of disaster victims. ? Mobilizing and positioning emergency supplies, equipment and personnel; including time-sensitive operations such as search and rescue, evacuation, emergency medical care, food and shelter programs, and bringing damaged services and systems back online. Recovery ? Actions that begin after the disaster, when urgent needs have been met. Recovery actions are designed to put the community back together ? Include repairs to roads, bridges, and other public facilities, restoration of power, water and other municipal services, and other activities to help restore normal operations to a community. Mitigation ? Activities that prevent a disaster, reduce the chance of a disaster happening, or lessen the damaging effects of unavoidable disasters and emergencies. ? Includes engineering solutions such as dams and levees; land-use planning to prevent development in hazardous areas; protecting structures through sound building practices and retrofitting; acquiring and relocating damaged structures; preserving the natural 9 environment to serve as a buffer against hazard impacts; and educating the public about hazards and ways to reduce risk. Emergency management process needs the cooperation of all individuals, groups, and communities to be successful. When a major disaster happens, emergency management agencies from all over the world work with governments and non-governmental organizations in an effort to decrease the impact of the disaster. Humanitarian organizations such as American Red Cross, CARE USA, Catholic Relief Services, International Committee of the Red Cross, International Federation of Red Cross and Red Crescent Societies, International Rescue Committee, UNICEF, World Bank, and World Food Program are among the organizations that work with different national organizations inside the affected countries to provide humanitarian aids. In the United States, the federal emergency management agency (FEMA) is the main agency to deal with emergencies. They work in partnership with other organizations that are part of the national emergency management system. These partners include state and local emergency management agencies, 27 other federal agencies and the American Red Cross. More details on FEMA?s structure and operations are introduced in the following section. 1.3 Federal Emergency Management Agency Federal emergency management agency (FEMA) is the main organization responsible for dealing with the federal level emergencies in the United States. It was initially created in 1979 as an independent organization but On March 1st, 2003 FEMA became part of the U.S. department of homeland security (DHS) 10 along with 22 other government agencies. FEMA is a relatively small agency with around 2,600 full time employees but it can mobilize nearly 7000 temporary disaster assistance employees to respond to disasters. Besides the headquarters in Washington D.C., FEMA has ten regional offices across the country to coordinate with its state and local government counterparts and with nonprofit and for-profit organizations. The primary mission of FEMA is ?To reduce the loss of life and property and protect the Nation from all hazards, including natural disasters, acts of terrorism, and other man-made disasters, by leading and supporting the Nation in a risk-based, comprehensive emergency management system of preparedness, protection, response, recovery, and mitigation.? (www.fema.gov) FEMA?s strategic plan for fiscal years 2008-2013 declares the vision of the organization as ?The nation?s preeminent emergency management and preparedness agency?. The plan establishes strategic goals, objectives, and strategies to fulfill FEMA?s vision. The strategic goals of the agency are to: 1. Lead an integrated approach that strengthens the nation?s ability to address disasters, emergencies, and terrorist events 2. Deliver easily accessible and coordinated assistance for all programs 3. Provide reliable information at the right time for all users 4. FEMA invests in people and people invest in FEMA to ensure mission success 5. Build public trust and confidence through performance and stewardship 11 One of the important documents that define the principles, roles, and structures of FEMA is the national response framework (NRF). NRF replaced its older version called the national response plan on March 22, 2008. NRF presents the guiding principles that enable all response partners to prepare for and provide a unified national response to disasters and emergencies. It describes how communities, tribes, states, the federal government, private-sectors, and nongovernmental partners work together to coordinate national response. Following the guidelines of NRF are essential to establish a comprehensive, national, all-hazards approach for disaster response in the United States. NRF main documents are supplemented by important annexes called emergency support functions (ESF). The ESFs provide the structure for coordinating federal interagency support for a federal response to an emergency. They are mechanisms for grouping functions most frequently used to provide federal support to states and federal-to-federal support, both for declared disasters and emergencies under the Stafford act and for non-Stafford act incidents. Table 1.1 gives a summary of the 15 ESFs currently present in the NRF. More information on the national response framework including documents, annexes, references and briefings/trainings can be accessed through the NRF resource center at www.fema.gov/nrf . 12 Table 1.1 Emergency Support Function Annexes of the National Response Framework Emergency Support Function Scope ESF 1 Transportation Aviation management and control; Transportation safety Restoration/recovery of transportation infrastructure; Movement restrictions; Damage and impact assessment ESF 2 Communications Coordination with telecommunications and information technology industries; Restoration and repair of telecommunications infrastructure. Protection, restoration, and sustainment of national cyber and information technology resources; Oversight of communications within the Federal incident management and response structures ESF 3 Public Works and Engineering Infrastructure protection and emergency repair; Infrastructure restoration; Engineering services and construction management; Emergency contracting support for life-saving and life-sustaining services ESF 4 Firefighting Coordination of Federal firefighting activities; Support to wild land, rural, and urban firefighting operations ESF 5 Emergency Management Coordination of incident management and response efforts; Issuance of mission assignments; Resource and human capital; Incident action planning; Financial management ESF 6 Housing, and Human Services Mass care; Emergency assistance; Disaster housing; Human services ESF 7 Logistics Management Comprehensive, national incident logistics planning, management, and sustainment capability; Resource support (facility space, office equipment and supplies, contracting services, etc.) ESF 8 Public Health Public health; Medical and Mental health services; Mass fatality management ESF 9 Search and Rescue Life-saving assistance Search and rescue operations ESF 10 Hazardous Materials Oil and hazardous materials (chemical, biological, radiological, etc.) response; Environmental short- and long-term cleanup ESF 11 Agriculture and Natural Resources Nutrition assistance; Animal and plant disease and pest response; Food safety and security; Natural and cultural resources and historic properties protection and restoration; Safety and well-being of household pets ESF 12 Energy Energy infrastructure assessment, repair, and restoration; Energy industry utilities coordination; Energy forecast ESF 13 Public Safety and Security Facility and resource security; Security planning and technical resource assistance; Public safety and security support; Support to access, traffic, and crowd control ESF 14 Long- Term Recovery Long-term community recovery assistance to States, local governments, and the private sector Analysis and review of mitigation program implementation ESF 15 External Affairs Emergency public information and protective action guidance; Media and community relations; Congressional and international affairs; Tribal and insular affairs 13 Emergency support function 7 is the logistics management and resource support annex that describes the roles and responsibilities of FEMA and general services administration (GSA) to jointly manage a supply chain that provides relief commodities to the victims. Based on ESF 7, FEMA is the primary agency for logistics management and is responsible for: ? Material management that includes determining requirements, sourcing, ordering and replenishment, storage, and issuing of supplies and equipment. ? Transportation management that includes equipment and procedures for moving material from storage facilities and vendors to incident victims, particularly with emphasis on the surge and sustainment portions of response. Transportation management also includes providing services to requests from other federal organizations. ? Facilities management that includes the location, selection, and acquisition of storage and distribution facilities. These facilities include logistics centers, mobilization centers, and federal operations staging areas. ? Personal property management and policy and procedures guidance for maintaining accountability of material and identification and reutilization of property acquired to support a federal response operation. ? Management of electronic data interchange to provide end-to-end visibility of response resources. 14 ? Planning and coordination with internal and external customers and other supply chain partners in the federal and private sectors for improving the delivery of goods and services to the customer. The next section introduces the components of FEMA?s logistics operations and describes the structure of FEMA?s supply chain. 1.4 FEMA?s Logistics Supply Chain FEMA has a complicated and special structure for its supply chain. There are seven main components in the FEMA?s supply chain to provide relief commodities for disaster victims that are briefly described here: 1. FEMA Logistics Centers (LC) - permanent facilities that receive, store, ship, and recover disaster commodities and equipment. FEMA has a total of 9 logistics centers: ? Four continental United States centers containing general commodities located at Atlanta, Georgia; Ft. Worth, Texas; Frederick, Maryland; and Moffett Field, California. ? Three off-shore centers containing general commodities located in Hawaii, Guam, and Puerto Rico. ? Two Continental United States centers containing special products such as computers, office electronic equipment, medical and pharmaceutical caches located in Cumberland, Maryland and Berryville, Virginia. 15 Examples of disaster relief commodities include ice, water, meals ready to eat (MREs), blankets, cots, flashlights, tarps, sleeping bags and tents. Disaster relief equipments include emergency generators, personal toilet kits, and refrigerated vans. 2. Commercial Storage Sites (CSS) - permanent facilities that are owned and operated by private industry and store commodities for FEMA. Freezer storage space for ice is an example. 3. Other Federal Agencies Sites (VEN) - representing vendors from whom commodities are purchased and managed. Examples are the defense logistics agency (DLA) and the general services administration (GSA). 4. Mobilization (MOB) Centers - temporary federal facilities in theater at which commodities, equipment and personnel can be received and pre- positioned for deployment as required. In MOBs commodities remain under the control of the FEMA logistics headquarter and can be deployed to multiple states. MOBs are generally projected to have the capacity to hold 3 days of supply commodities. 5. Federal Operational Staging Areas (FOSAs) - temporary facilities at which commodities, equipment and personnel are received and pre-positioned for deployment within one designated state as required. Commodities are under the control of the operations section of the joint field office (JFO) or the regional response coordination center (RRCC). Commodities are usually 16 being supplied from MOB centers, logistics centers or direct shipments from vendors. FOSAs are generally projected to hold 1 to 2 days of commodities. 6. State Staging Areas (SSA) - temporary facilities in the affected state at which commodities, equipment and personnel are received and pre-positioned for deployment within that state. Title transfers for delivered federal commodities and cost sharing are initiated in SSAs. 7. Points of Distribution (PODs) Sites - temporary local facilities in the disaster area at which commodities are distributed directly to disaster victims. PODs are operated by the affected state. Figure 1.5 better illustrates this structure. At the top of the pyramid there are 3 types of facilities namely FEMA logistics centers, commercial storage sites, and other federal agencies or vendors. These permanent facilities store and ship commodities and equipment and are considered as the ?sources? of the chain. The mobilization centers, the federal operational staging areas, and the state staging areas are 3 types of facilities that mainly play the role of transshipment points. These are temporary facilities at which commodities, equipment and personnel are received and pre-positioned for deployment to the lower levels. At the end, point of distribution sites are temporary local facilities at which commodities are received and distributed directly to the disaster victims. The PODs can be local schools, churches, or even large parking lots in the affected area. 17 Figure 1.5 FEMA?s Supply Chain Structure Even this simplified presentation of the FEMA?s logistics supply chain indicates the complex structure of the system. Finding the optimal sites for 4 levels of temporary facilities is a complicated location finding problem. Delivering several types of relief commodities to the disaster victims is a multicommodity capacitated network flow problem. Optimizing the movement of vehicles in the network is a dynamic vehicle routing problem with mixed pickup and delivery operations. Usually more than one transportation mode is used in disaster response operations which makes the problem a multimodal transportation problem. Other characteristics that make the problem unique include, but are not limited to, importance of quick response and fast delivery, Logistics Centers MOB Centers State Staging Area POD Commercial Storage Sites State Staging Area POD POD POD Vendors FOSAs 18 shortage of supply versus overwhelming demands, insufficient capacity of the facilities and transportation system, and dynamic environment of the emergency situations. 1.5 Motivation and Objective of the Research In today?s society that disasters seem to be striking all corners of the United States and the globe, the importance of the emergency management is undeniable. Much human loss and unnecessary destruction of infrastructure can be avoided with more foresight and specific planning as well as a precise execution. In a world where resources are stretched to the limit and the question of humanitarian relief seems too often to be tied with economical considerations, better designs and operations are urgently needed to help save thousands of lives and millions of dollars. The question is how to respond to natural disasters in the most efficient manner to minimize the loss of life and maximize the efficiency of the rescue operations. In case of these emergencies various organizations often face significant problems of transporting large amounts of many different commodities including food, clothing, medicine, medical supplies, machinery, and personnel from different points of origin to different destinations in the disaster areas. The transportation of supplies and relief personnel must be done quickly and efficiently to maximize the survival rate of the affected population and minimize the cost of such operations. Federal emergency management agency (FEMA) is the primary organization for preparedness and response to the federal level disasters in the 19 United States. Unfortunately, inadequate response to hurricanes Katrina and Rita showed the critical need for better mechanisms in emergency operations. Initial research in this area shows that this is an emerging field and there are great potentials for research in emergency logistics and disaster response. FEMA has a very complex logistics structure to provide the disaster victims with critical items after a disaster strike which involves multiple organizations and spreads across the country. The goal of this research is to develop a comprehensive model that describes the integrated logistics operations in response to natural disasters at the operational level. The proposed mathematical model integrates three main components. First, it controls the flow of several relief commodities from sources through the supply chain until they are delivered to the hands of the recipients. Second, it considers a large-scale unconventional vehicle routing problem with mixed pickup and delivery schedules for multiple transportation modes. And third, following the FEMA?s complex logistics structure, a special facility location problem is considered that involves four layers of temporary facilities at the federal and state levels. Such integrated model provides the opportunity for a centralized operation plan that can effectively eliminate the delays and assign the limited resources in a way that is optimal for the entire system. 1.6 Contributions of the research Emergency response is a dynamic and very time sensitive operation. This research offers an integrated model that not only considers details such as multimodal vehicle routing and pick up or delivery schedules; but also considers 20 finding the optimal locations for the temporary facilities as well as considering the capacity constraint for each facility and the transportation system. A mathematical model at the operational level is presented that can be used in the critical hours and days immediately after the disaster strikes. Such a model is a unique tool that can also be used at strategic level or planning level analysis. It is a very complicated task and up to date, there is no study in the literature that has addressed this problem sufficiently. This research also aims at developing optimization algorithms and heuristics to solve the proposed model and find applicable solutions to decrease human sufferings in the most economically sensible way. The algorithms need to be fast so that the results can be used in the initial response phase and also as the situation changes in the highly dynamic environment after the disaster. Also, in this research a comprehensive set of numerical case studies and sensitivity analysis are performed. In-depth analyses of different aspects of the proposed mathematical model are provided in order to better illustrate the capabilities of the model and also examine model?s sensitivity in various circumstances. These analyses are intended to introduce a general framework for researchers and practitioners. The findings of these analyses may or may not be directly applicable for other specific disaster response scenarios; however, they provide a general study framework for modeling and analysis of each specific disaster scenario that can be adopted by other researchers and practitioners. In other words, this research extends the state-of-the-art by presenting an integrated model at the operational level that describes the details of the supply 21 chain operations in major emergency management agencies such as FEMA, in response to immediate aftermath of a large scale disaster. Development of fast and efficient solution algorithms and heuristics for the proposed model is the other major contribution of this research. The simulations and sensitivity analysis provided in this research can be used as a framework to follow by other researchers and practitioners. 1.7 Organization of the dissertation After this introduction, previous works in the fields of logistics and disaster relief operations are reviewed in chapter 2. The specific problem to be dealt with in this research is introduced in chapter 3 and then the mathematical formulation of the model is presented. Chapter 4 offers a set of preliminary numerical examples to evaluate the model and help better understand the mechanics of the model. In chapter 5, solution approaches are summarized and a number of heuristic solution algorithms are proposed to solve the different parts of the proposed model. Chapter 6 is dedicated to simulation and in-depth analysis of the proposed model and sensitivity analysis of its parameters. In chapter 7, the prepositioning of relief supplies and equity constraints are discussed. Finally in chapter 8, a summary of this dissertation is presented and some suggestions for future research are discussed. 22 Chapter 2: Literature review In this chapter, first in section 2.1 some definitions of supply chain and supply chain management (SCM) in commercial sector are introduced then some of the researches that reviewed the supply chain studies are summarized. Then a brief introduction to facility location problem is presented in section 2.2 and vehicle routing problem in section 2.3 are presented as two main elements of the supply chain and logistics modeling. In section 2.4, the similarities and differences between commercial supply chain and disaster response logistics are reviewed. In section 2.5, some studies specific to modeling and optimization of logistics in disaster response are provided. Finally in section 2.6, a summary of previews works in this area is presented with the emphasis on the gaps in the literature that needs to be filled. 2.1 Supply Chain Management Definition of supply chain management differs across authors from different fields and there is no explicit and universal description of SCM or its activities in the literature (Tan 2001). The literature is full of buzzwords such as integrated purchasing strategy, integrated logistics, supplier integration, buyer- supplier partnerships, supply base management, strategic supplier alliances, supply chain synchronization and supply chain management, to address elements or stages of this phenomenon (New, 1997; La Londe and Masters, 1994). For example Harland (1996) described supply chain management as managing business activities and relationships (1) internally within an 23 organization, (2) with immediate suppliers, (3) with first and second-tier suppliers and customers along the supply chain, and (4) with the entire supply chain. Scott and Westbrook (1991) and New and Payne (1995) described supply chain management as the chain linking each element of the manufacturing and supply process from raw materials through to the end user, including several organizational boundaries. SCM begins with the extraction of raw materials or minerals from the earth, through the manufacturers, wholesalers, retailers, and the final users. Where appropriate, supply chain management also includes recycling or re-use of the products or materials. Another definition of supply chain management emerges from the transportation and logistics literature of the wholesale and retail industry, emphasizing the importance of physical distribution and integrated logistics. There is no doubt that logistics is an important function of business and is evolving into strategic supply chain management (New and Payne, 1995). In this definition, the physical transformation of the products is not a critical component of supply chain management. Its primary focus is the efficient physical distribution of final products from the manufacturers to the end users in an attempt to replace inventories with information and reduce transportation costs. The definition of supply chain seems to be more common across authors than the definition of supply chain management (Mentzer et al. 2001). La Londe and Masters (1994) proposed that the supply chain is a set of firms that pass materials forward. Eksioglu (2002) defined the supply chain as an integrated process where different business entities such as suppliers, manufacturers, 24 distributors, and retailers work together to plan, coordinate, and control the flow of materials, parts, and finished goods from suppliers to customers. Several independent firms can be involved in manufacturing a product and placing it in the hands of the end user in a supply chain. For example raw material and component producers, product assemblers, wholesalers, retailer merchants and transportation companies are all members of the supply chain. Beamon (1998) defined supply chain as an integrated manufacturing process where raw materials are converted into final products, then delivered to customers. At its highest level, a supply chain is comprised of two basic integrated processes: (1) the production planning and inventory control process, and (2) the distribution and logistics process. These processes define the basic framework for the conversion and movement of raw materials into final products. Figure 2.1 illustrates a simplified picture of the supply chain process. Figure 2.1 Supply chain process (adopted from Beamon 1998) The production planning and inventory control process includes the manufacturing and storage sub-processes and their interfaces. More specifically, production planning describes the design and management of the entire manufacturing process including raw material scheduling and acquisition, Suppliers Manufacturers Storage Facilities Transportation Distribution Centers Retailers Production planning and Inventory Control Distribution and Logistics 25 manufacturing process design and scheduling, and material handling design and control. Inventory control describes the design and management of the storage policies and procedures for raw materials, work-in-process inventories, and usually, final products. The distribution and logistics process determines how products are retrieved and transported from the storage warehouse to retailers. These products may be transported to retailers directly, or may be shipped to distribution facilities first and then being delivered to the retailers. This process includes the management of inventory retrieval, transportation, and final product delivery. These processes interact with one another to produce an integrated supply chain. The design and management of these processes determine the extent to which the supply chain works as a unit to meet required performance objectives. Usually in commercial supply chain, the objective is to minimize cost. However, some have considered a combination of cost and customer service as the objective of the commercial supply chains. For many years, researchers and practitioners have concentrated on the individual processes and entities within the supply chains. However, the recent trend is to model and optimize SC as a single unified entity. In this approach, operations research (OR) techniques have shown to be a very useful tool among researchers and practitioners. Typically, a SC model tries to determine ? the transportation modes to be used, ? the suppliers to be selected, ? the amount of inventory to be held at various locations in the chain, 26 ? the number of warehouses to be used, and ? the location and capacities of these warehouses For a more comprehensive review of models and methods in supply chain design and analysis, readers are referred to Beamon (1998) and Tan (2001). In the following sections, some of the elements of SCM that can be applied in disaster response logistics are introduced in more details. 2.2 Facility Location Problem One of the most important problems in supply chain management is deciding where to locate new facilities such as factories, warehouses, distribution centers or retailers to support the material flow through an efficient distribution system. The general facility location problem can be stated as: for a given set of facility locations and a set of customers who are served from these facilities, find: ? Which facilities should be used ? Which customers should be served from which facilities so as to minimize the total cost of serving all the customers The development and acquisition of a new facility is typically a costly and time-consuming project. Before a facility can be purchased or constructed, good locations must be identified, appropriate facility capacity specifications must be determined, and large amounts of capital must be allocated. While the objectives driving a facility location decision depend on the firm or government agency, the 27 high costs associated with this process make almost any location project a long- term investment. A vast literature has developed out of the broadly based interest in facility location problem over the last four decades (Daskin 1995, Drezner and Hamacher 2002). Operations research practitioners have developed a number of mathematical programming models to represent a wide range of location problems. Several different objective functions have been formulated to consider numerous applications. Unfortunately, the resulting models can be extremely difficult to solve to optimality (most problems are classified as NP-hard); many of the problems require integer programming formulations. The p-median problem, covering problem, and p-center problem are three classic forms of facility location problem that are introduced in the following subsections. For a comprehensive bibliography of more recent studies in discrete location finding problem refer to ReVelle et al (2008). 2.2.1 P-Median Problem One important way to measure the effectiveness of a facility?s location is by determining the average distance traveled by those who visit it. As average travel distance increases, facility accessibility decreases, and thus the location's effectiveness decrease. An equivalent way to measure location effectiveness when demands are not sensitive to the level of service is to weight the distance between demand nodes and facilities by the associated demand quantity and calculate the total weighted travel distance between demands and facilities. Then, the problem is to selects the best p sites among a range of possible locations with the objective 28 of minimizing total demand-weighted travel distance between demand nodes and selected facilities. The key decisions are where to locate the p facilities and which facility should serve each demand node. The inputs are the demands (or weights) iw at each node Ii ? , the distances dij between each demand node Ii ? and each candidate facility site Jj ? and p, the maximum number of facilities to be located. The mathematical formulation of p-median problem is as follow: xj = 1 if a facility is located at candidate node Jj ? and 0 otherwise yij = 1 if demand node Ii ? is assigned to facility at candidate node Jj ? 0 otherwise. Minimize ?? ? ?Jj Ii ijiji ydw (2.1) Subject to Iiy Jj ij ??=? ? 1 (2.2) JjIixy jij ?????? ,0 (2.3) Jjpx Jj ij ???? ? (2.4) { } JjIiyx ijj ????? ,1,0, (2.5) The objective function (2.1) minimizes the demand-weighted total distance. Since the demands are known and the total demand is fixed, this is equivalent to minimizing the demand-weighted average distance. Constraints (2.2) ensure that each demand node is assigned, while constraints (2.3) stipulate 29 that the assignments can only be made to open facilities. Constraint (2.4) states that a maximum of p facilities are to be opened. Constraints (2.5) are standard integrality constraints. 2.2.2 Covering Problem The p-median problem described above can be used to locate a wide range of public and private facilities. For some facilities, however, selecting locations which minimize the average distance traveled may not be appropriate. Suppose, for example, that a city is locating emergency service facilities such as fire stations or ambulances. The critical nature of demands for service will dictate a maximum ?acceptable? travel distance or time. Such facilities will thus require a different measure of location efficiency. To locate such facilities, the key issue is the ?coverage?. A demand is said to be covered if it can be served within a specified time. The literature on covering problem is divided into two major segments, that in which coverage is required and that in which it is optimized. Two covering problems which illustrate the distinction are the location set covering problem and the maximal covering problem. We introduce both problem classes. For a more complete review of covering problems refer to Schilling et al (1993). In the set covering problem, the objective is to minimize the cost of facility location such that a specified level of coverage is obtained. The mathematical formulation of set covering problem is as follow: cj = fixed cost of locating a facility at node j S = maximum acceptable distance or travel time 30 Ni = set of facility sites j within acceptable distance of node i ( { }SdjN iji ?= ) Xj = 1 if a facility is located at candidate node Jj ? and 0 otherwise Minimize ? ?Jj jj Xc (2.6) Subject to iX iNj j ??? ? 1 (2.7) { } jX j ?? 1,0 (2.8) The objective function (2.6) minimizes the cost of facility location. In many cases, the costs cj are assumed to be equal for all potential facility sites j, implying an objective equivalent to minimizing the number of facilities located. Constraint (2.7) requires that all demand nodes i have at least one facility located within the acceptable service distance. Note that this formulation makes no distinction between nodes based on demand size. Each node, whether it contains a single customer or a large portion of the total demand, must be covered regardless of the cost. If the coverage distance S is small, relative to the spacing of demand nodes, the coverage restriction can lead to a large number of facilities being located. Additionally, if an outlying node has a small demand, the cost/benefit ratio of covering that demand can be extremely high. In many practical applications, decision makers find that their allocated resources are not sufficient to build the facilities dictated by the desired level of coverage. In other words, the goal of coverage within distance S may be infeasible with respect to construction resources. In such cases, location goals must be shifted so that the available resources are used to give the desired level of 31 coverage to as many customers as possible. This new objective is that of the maximal covering problem. Specifically, the maximal covering problem seeks to maximize the amount of demand covered within the acceptable service distance S by locating a fixed number of facilities: Xj = 1 if a facility is located at candidate node Jj ? and 0 otherwise Zi = 1 if a demand at node Ii ? is covered and 0 otherwise Minimize ? i ii Zh (2.9) Subject to iXZ iNj ji ?? ? ? (2.10) ipX j j ??? (2.11) { } jiZX ij ,1,0, ?? (2.12) The objective (2.9) is to maximize the amount of demand covered. Constraint (2.10) determines which demand nodes are covered within the acceptable service distance. Each node i can only be considered covered (with Zi = 1) if there is a facility located at some site j which is within S of node i (i.e., if Xj = 1 for some iNj ? ). If no such facility is located, the right hand side of constraint (2.10) will be zero, thus forcing Zi to be zero. Constraint (2.11) limits the number of facilities to be located to a fixed number p. 32 2.2.3 P-Center Problem Another problem class which avoids the set covering problem's potential infeasibility is the class of p-center problems. In such problems, we require coverage of all demands, but we seek to locate a given number of facilities in such a way that minimizes coverage distance. Rather than taking an input coverage distance S, this model determines endogenously the minimal coverage distance associated with locating p facilities. The p-center problem is also known as the minimax problem, as we seek to minimize the maximum distance between any demand and its nearest facility. If facility locations are restricted to the nodes of the network, the problem is a vertex center problem. Center problems which allow facilities to be located anywhere on the network are absolute center problems. The following additional decision variable is needed in order to formulate the p-center problem: D = maximum distance between a demand node and the nearest facility. The resulting integer programming formulation of the P-center problem is as following: Minimize D (2.13) Subject to ipX j j ??? (2.14) iY j ij ?=? 1 (2.15) jiXY jij ,0 ??? (2.16) 33 iDYd j ijij ??? (2.17) { } jiYX ijj ,1,0, ?? (2.18) The objective function (2.13) is simply to minimize the maximum distance between any demand node and its nearest facility. Constraints (2.14) limits the maximum number of open facilities to p. constraints (2.15) enforces each demand point to be assigned to a facility and constraints (2.16) make sure that demands are assigned only to selected facilities. Constraint (2.17) defines the maximum distance between any demand node i and the nearest facility j. Finally, constraints (2.18) are integrality constraints for the decision variables. In addition to three classes introduced here, several alternate formulations of the facility location problem are proposed by researchers over the years. For a bibliography of recent studies refer to ReVelle et al. (2008). 2.3 Vehicle Routing Problem The vehicle routing problem (VRP) is a generic name given to a whole class of problems in which a set of routes for a fleet of vehicles based at one or several depots must be determined for a number of geographically dispersed cities or customers. The VRP arises naturally as a central problem in the fields of transportation, distribution and logistics. Usually, the objective of the VRP is to deliver a set of customers with known demands on minimum-cost vehicle routes originating and terminating at a depot. In some market sectors, transportation means a high percentage of the value added to goods. Therefore, the utilization of modeling and optimization methods for transportation often results in significant 34 savings ranging from 5% to 20% in the total costs, as reported in Toth and Vigo (2002). The VRP is a well known integer programming problem which falls into the category of NP-Hard problems, meaning that the computational effort required for solving this problem increases exponentially with the problem size. This difficult combinatorial problem conceptually lies at the intersection of these two well-studied NP-Hard problems: ? The Traveling Salesman Problem (TSP): If the capacity of the vehicles is infinite, we can get an instance of the multiple traveling salesman problem (MTSP). An MTSP instance can be transformed into an equivalent TSP instance by adding to the graph k-1 (k being the number of routes) additional copies of node 0 and its incident edges. ? The Bin Packing Problem (BPP): The question of whether there exists a feasible solution for a given instance of the VRP is an instance of the BPP. The decision version of this problem is conceptually equivalent to a VRP model in which all edge costs are taken to be zero (so that all feasible solutions have the same cost). Three basic approaches have been proposed for modeling VRP in the literature (Toth and Vigo 2002). The models of the first type, known as vehicle flow formulation, use binary integer variables associated with each arc of the network, which shows if an specific arc is traversed by a vehicle or not. These models are often used for basic versions of VRP. They are particularly useful for cases in which the cost of the solution can be expressed as the sum of the costs 35 associated with the arcs. On the other hand, vehicle flow models cannot be used to deal with many practical issues; for instance, when the cost of a solution depends on the sequence of traversed arcs or when the cost depends on the type of vehicle that is assigned to a route. The second approach to VRP modeling is called commodity flow formulation. In this type of model, additional integer variables are associated with arcs that represent the flow of the commodities along the paths traveled by the vehicles. In some recent studies, these models have been used as a basis to solve for the exact solutions of capacitated VRP. In the third approach to VRP modeling, the decision variables are the feasible routes for the vehicles. These models produce an exponential number of binary variables each associated with a feasible route. Then the VRP is formulated as a set partitioning problem that tries to select a set of routes with minimum cost which serves each costumer once and also satisfies the additional constraints. Main advantage of this type of model is that it allows for extremely general route costs. For example, route costs can be nonlinear or can depend on the vehicle type or sequence of nodes visited. Also, the linear relaxation of these models usually provides a tighter bound than the previous models. However, these models usually require enumerating the feasible routes which needs extensive preprocessing and results in a very large number of variables. 2.3.1 VRP Mathematical Formulation As mentioned above, vehicle flow based formulation is one of the approaches to model the VRP. Following formulation is an example for the base 36 case of uncapacitated multi-vehicle single depot vehicle routing problem. The decision variables vijx which are binary indicate whether vehicle v travels from point i to point j, vijx =1, or not vijx =0 Minimize ??? i j v v ijij xc (2.19) Subject to jx v i v ij ?=?? 1 (2.20) ix v j v ij ?=?? 1 (2.21) vNpxx j v pj i v ip ???=? ?? ,0 (2.22) vx j v j ??? 10 (2.23) ( ) vjixvij ,,1,0 ?? (2.24) SX ? (2.25) The objective is to minimize the total travel cost (or distance) by all vehicles. Constraints (2.20) through (2.22) require that only one vehicle enters each node and that the same vehicle exits that node. Constraints (2.23) insure that each vehicle leaves the depot only once. The last condition which is imposed on the matrix X prohibits sub-tours that do not contain the depot. There are several possible ways to fulfill this condition, for example S might be composed of sub- tour breaking constraints for each vehicle. S can be defined as the union of sets Sv defined by: 37 ? ? ? ? ? ? ??= ?? ? ?Qi Qj v ij v ijv QsubsetnonemptyallforQxxS 1: (2.26) If each customer has a demand of di units and each vehicle has a capacity of Kv, then the capacitated VRP can be formulated by adding the following capacity constraints to the base formulation: vKxd v i j v iji ????? ? ??? ?? ? (2.27) 2.3.2 VRP Variants Usually, real world vehicle routing problems are much more sophisticated than the base case VRP introduced above. Over the years, researchers have proposed variants of VRP by adding some constraints to the base case VRP formulation. Here, a list of well-known VRP variants is summarized: ? Capacitated VRP (CVRP): Every vehicle has a limited capacity ? Distance-Constrained VRP (DCVRP): The maximum tour length is limited ? Multiple Depot VRP (MDVRP): The vendor uses many depots to supply the customers ? VRP with Pick-Up and Delivering (VRPPD): Customers may return some goods to the depot or other customers ? Split Delivery VRP (SDVRP): The customers may be served by different vehicles ? VRP with time windows (VRPTW): Every customer has to be supplied within a certain time window 38 ? Periodic VRP (PVRP): The deliveries may be done in some consecutive days ? Stochastic VRP (SVRP): Some values, such as number of customers, their demands, service time or travel time, are stochastic variables There are several survey papers on the VRP, VRP variants, and their solution algorithms and techniques. A classification of the problem was given in Desrochers et al.(1990). Laporte and Nobert (1987) presented a survey of exact methods to solve VRP. Other surveys that provided exact and heuristic methods were presented by Christofides, Mingozzi, and Toth (1981), Magnanti (1981), Bodin et al.(1983), Fisher (1994), Laporte (1992), Toth and Vigo (2002). An annotated bibliography was proposed by Laporte (1997). A book on the subject was edited by Golden and Assad (1988). 2.4 Commercial Supply Chain vs. Emergency Response Logistics Immediately after a disaster, humanitarian organizations often face significant problems of transporting large amounts of many different commodities including food, clothing, medicine, medical supplies, machinery, and personnel from several origins to several destinations inside the disaster area. The transportation of supplies and relief personnel must be done quickly and efficiently to maximize the survival rate of the affected population and minimize the cost of such operations. When it comes to efficiency of supply deliveries, the modeling and optimization techniques established in commercial supply chain management seem to be the most relevant approach. For instance, some of the quickest emergency assistance to the victims of hurricane Katrina did not come from the 39 American Red Cross or FEMA, it came from Wal-Mart. Millions of affected or displaced people waited for days as agencies struggled to provide assistance. Wal- Mart moved faster than traditional emergency aid groups mainly because the retail giant had mastered the fundamentals of logistics and supply chain management (Dimitruk 2005). More recently, some studies such as (Beamon 2004; Thomas and Kopczak, 2005; Van Wassenhove, 2006; Oloruntoba and Gray, 2006; Thomas, 2007), emphasized that some supply chain concepts share similarities to emergency logistics and therefore some tools and methods developed for commercial supply chains can be successfully adapted in emergency response logistics. Using commercial supply chain techniques in disaster management is still in its infancy. Beamon (2004) and Thomas (2005) have compared the current state of supply chain management capabilities within humanitarian organizations with that of the commercial sector in the 1970s and 1980s. At that time, the commercial sector just began to realize the strategic advantages and significant improvements supply chain management could offer in effectiveness and efficiency. This led to extensive research in the area of supply chain and logistical analysis but those quantitative methods and principles are rarely applied to humanitarian operations on the verge of disasters. The partial reason is the difference in the strategic goals of commercial supply chain with goals of disaster response logistics. The main goal in commercial supply chain is to minimize the cost or maximize the profit of 40 operations. Actions are justified if they increase the profit but are not perused if their cost is more than their profit. However, humanitarian organizations are mostly non-profit organizations with the idea of providing critical services to the public in order to minimize the pain and sufferings, for example after a natural disaster. One major difference between the two types of chains is the demand pattern. For many commercial supply chains, the external demand for products is comparatively stable and predictable. Often, for the commercial chain, the demands seen from warehouses occur from established locations in relatively regular intervals. However, the demands in the relief chain are emergency items, equipment, and personnel. More importantly, those demands occur in irregular amounts and at irregular intervals and occur suddenly, such that the locations are often completely unknown until the demand occurs. Beamon (2004) suggests other specific characteristics of disaster response logistics that differentiate them from traditional commercial supply chains. These include: ? Zero lead-time that dramatically affects inventory availability, procurement, and distribution ? High stakes (often life-and-death) that requires speed and efficiency ? Unreliable, incomplete, or non-existent supply and transportation infrastructure ? Many relief operations are naturally ad hoc, without effective monitoring and control 41 ? Variable levels of technology is available depending on the disaster area Table 2.1 compares some of the differences between commercial and humanitarian supply chains. Table 2.1- Commercial supply chains vs. humanitarian relief chains (Beamon 2004) Characteristic Commercial Chain Humanitarian Relief Chain Strategic Goals Typically to produce high quality products at low cost to maximize profitability Minimize loss of life and alleviate suffering Distribution Network Configuration Well-defined methods for determining the number and locations of distribution centers. Challenging due to the nature of the unknowns (locations, type and size of events, politics, and culture) Demand Type Commercial Products Emergency Supplies, equipment and Personnel Lead Time Lead time determined by the supplier-manufacturer-DC-retailer Zero time between the occurrence of the demand and the need for the demand Inventory Control Utilizes well-defined methods for determining inventory levels based on lead time, demand and target customer service levels Inventory control is challenging due to the high variations in lead times, demands, and demand locations Information System Generally well-defined, using advanced technology Information is often unreliable, incomplete or non-existent It is concluded that some of the concepts associated with commercial supply chains are directly applicable to humanitarian relief chains. However, future work must develop methods that specifically address the challenges 42 presented by characteristics unique to humanitarian relief and logistics of disaster response. 2.5 Logistics in Disaster Response Altay and Green (2006) surveyed the existing literature of emergency disaster management. They concluded that most of the disaster management research was related to social sciences and humanities literature. Refer to Hughes (1991) and http://www.geo.umass.edu/courses/geo510/index.htm for a comprehensive bibliography. That type of research focuses on subjects such as disaster results, sociological impacts on communities, psychological effects on survivors or rescue teams, and organizational design and communication problems. They observed that the existing literature is relatively light on disaster management articles that used operations research or management science (OR/MS) techniques to deal with the problem. However, they realized the literature trend that more studies are focusing on OR/MS techniques in recent years and emphasized the need for more research in future. In the following subsections, a summary of studies is presented that use OR/MS techniques to model and optimize the emergency disaster management activities. This is not an exclusive list of publication in the field and is only intended to focus on key studies in the past that successfully used techniques that are relevant to the subject of this dissertation. 43 2.5.1 Early Ages A number of authors have recognized the problem of emergency response management in its early ages. Kemball-Cook and Stephenson (1984) addressed the need for logistics management in relief operations for the increasing refugee population in Somalia. Ardekani and Hobeika (1988) addressed the need of logistics management in relief operations for the 1985 Mexico City earthquake. Knott (1987) developed a linear programming model for the bulk food transportation problem and the efficient use of the truck fleet to minimize the transportation cost or to maximize the amount of food delivered (single commodity, single modal network flow problem). In another article, Knott (1988) developed a linear programming model using expert knowledge for the vehicle scheduling of bulk relief of food to a disaster area. Ray (1987) developed a single-commodity, multi-modal network flow model on a capacitated network over a multi-period planning horizon to minimize the sum of all costs incurred during the transport and storage of food aid. Brown and Vassiliou (1993) developed a real-time decision support system which uses optimization methods, simulation, and the decision maker?s judgment for operational assignment of units to tasks and for tactical allocation of units to task requirements in repairing major damage to public works following a disaster. The literature in the multi-commodity, multi-modal network flow problem was relatively sparse. Crainic and Rousseau (1986) developed an optimization algorithm based on decomposition and column generation principles to minimize the total operating and delay cost for multi-commodity, multi-modal freight 44 transportation when a single organization controls both the service network and the transportation of goods. Guelat et al. (1990) presented a multi-commodity, multi-modal network assignment model for the purpose of strategic planning to predict multi-commodity flows over a multi-modal network. The objective function to be minimized was the sum of total routing cost and total transfer costs. 2.5.2 Recent Studies Technology advancement in recent years has opened new doors for researchers. Haghani and Oh (1996) proposed a formulation and solution of a multi-commodity, multi-modal network flow model for disaster relief operations. Their model could determine detailed routing and scheduling plans for multiple transportation modes carrying various relief commodities from multiple supply points to demand points in the disaster area. They formulated the multi-depot mixed pickup and delivery vehicle routing problem with time windows as a special network flow problem over a time-space network. The objective was minimizing the sum of the vehicular flow costs, commodity flow costs, supply/demand storage costs and inter-modal transfer costs over all time periods. They developed two heuristic solution algorithms; the first was a Lagrangian relaxation approach, and the second was an iterative fix-and-run process. Their work is one of the few studies that can be implemented at the operational level. Barbarosoglu et al. (2002) focused on tactical and operational scheduling of helicopter activities in a disaster relief operation. They proposed a bi-level modeling framework to address the crew assignment, routing and transportation issues during the initial response phase of disaster management in a static manner. 45 The top level mainly involves tactical decisions of determining the helicopter fleet, pilot assignments and the total number of tours to be performed by each helicopter without explicitly considering the detailed routing of the helicopters among disaster nodes. The base level addresses operational decisions such as the vehicle routing of helicopters from the operation base to disaster points in the emergency area, the load/unload, delivery, transshipment and rescue plans of each helicopter in each tour, and the re-fueling schedule of each helicopter given the solution of the top level. Barbarosoglu and Arda (2004) developed a two-stage stochastic programming model for transportation planning in disaster response. Their study expanded on the deterministic multi-commodity, multi-modal network flow problem of Haghani and Oh (1996) by including uncertainties in supply, route capacities, and demand requirements. The authors designed 8 earthquake scenarios to test their approach on real-world problem instances. It is a planning model that does not deal with the important details that might be required at strategic or operational level. It does not address facility location problem or vehicle routing problem. Ozdamar et al. (2004) addressed an emergency logistics problem for distributing multiple commodities from a number of supply centers to distribution centers near the affected areas. They formulated a multi-period multi-commodity network flow model to determine pickup and delivery schedules for vehicles as well as the quantities of loads delivered on these routes, with the objective of minimizing the amount of unsatisfied demand over time. The structure of the 46 proposed formulation enabled them to regenerate plans based on changing demand, supply quantities, and fleet size. They developed an iterative Lagrangian relaxation algorithm and a greedy heuristic to solve the problem. Yi and Ozdamar (2007) proposed a model that integrated the supply delivery with evacuation of wounded people in disaster response activities. They considered establishment of temporary emergency facilities in disaster area to serve the medical needs of victims immediately after disaster. They used the capacity of vehicles to move wounded people as well as relief commodities. Their model considered vehicle routing problem in conjunction with facility location problem. The proposed model is a mixed integer multi-commodity network flow model that treats vehicles as integer commodity flows rather than binary variables. That resulted in a more compact formulation but post processing was needed to extract detailed vehicle routing and pick up or delivery schedule. They reported that post processing algorithm was pseudo-polynomial in terms of the number of vehicles utilized. In a recent study, Balcik and Beamon (2008) proposed a model to determine the number and locations of distribution centers to be uses in relief operations. They formulated the location finding problem as a variant of maximum covering problem when the demand estimations are available for a set of likely scenarios. Their objective function maximizes the total expected demand covered by the established distribution centers. They also solve for the amount of relief supplies to be stocked at each distribution center to meet the demands. Their study is one of the first to solve location finding problem in relief operation; 47 however, they do not consider the location problem as part of a supply chain network. Consequently, they cannot consider the interactions between optimal transportation of relief items from sources to the demand points and problem of finding optimal locations for distribution facilities. 2.6 Conclusions There are not many publications that directly applied network modeling and optimization techniques in disaster response. Among those studies, there is no model that has integrated the interrelated problems of large scale multicommodity multimodal network flow problem, vehicle routing problem with split mixed pickup and delivery, and optimal location finding problem with multiple layers. Also to the best of our knowledge, there is no mathematical model that describes the special structure of FEMA?s supply chain system. It is intended to fill some of these gaps in the following chapters of this research. After providing a more formal description of the problem, a mathematical model is proposed that considers the specific characteristics of the described problem. The proposed mathematical model is a comprehensive system that integrates all the above mentioned properties. Offering this large-scale mathematical formulation is a unique theoretical contribution by itself. Nevertheless, solving this large-scale integrated formulation for real-world size problems requires special considerations. This problem belongs to the NP-Hard class that is proven to be extremely time-consuming as the problem size grows. Offering fast and efficient solution algorithms and heuristics is another gap that is being addressed in this research. 48 Finding fast solution algorithms is especially important because it enables the real-time optimization and implementation of the proposed model. Extensive numerical and sensitivity analysis are required to evaluate the different aspects of the proposed model and solution algorithms. Through numerical case studies and simulation scenarios, it will be possible to fully test the performance of the model and the solution algorithms. The other important outcome of extensive numerical and sensitivity analysis will be the development of a set of general guidelines for practitioners, in order to model and solve similar case studies for their specific applications. 49 Chapter 3: Problem Description and Formulation In this chapter, first a complete description of the problem and its properties are provided in section 3.1. In section 3.2, the concept of time-space network is introduced which is very important in modeling the dynamic behavior of the described problem. Then the research approach to build a mathematical model for the problem is described in section 3.3 followed by the list of assumptions made in order to properly model the problem. In section 3.5, the details of mathematical formulation for the proposed model are presented. The notations, parameters and variables are defined in sections 3.5.1 through 3.5.3. The objective function of the optimization problem is formulated in section 3.5.4, followed by the formulation and description of the constraints of the problem in section 3.5.5. Finally, in section 3.6, a short form of the mathematical formulation is presented for the summary. 3.1 Problem Description The goal of the mathematical model is to orchestrate all the logistical components and tasks in the emergency response operations after a large scale disaster, in order to minimize the loss of life or human sufferings by rapid and efficient delivery of critical relief items to the victims in the disaster areas. Logistics planning in emergencies involves sending multiple relief commodities (e.g., medicine, water, food, equipment, etc) from a number of sources to several distribution points in the affected areas through a chain structure with some intermediate transfer nodes. The supplies may not be 50 available immediately but arrive over time. It is a difficult task to decide on the right type and quantity of relief items, the sources and destinations of the commodities, and also how to dispatch relief items to the recipients in order to minimize the pain and sufferings for the disaster victims. It is necessary to have a quick estimation of the demands during the initial response time. It is essential to know the types of required commodities, the amounts of each commodity per person or household, an estimation of the number of victims, and the geographical locations of the demands. The list of commodities includes but is not limited to water, food, shelter, electric generators, medical supplies, cots, blankets, tarps and clothing. Some of the demand items are one-time demands while others are recurring (e.g. tent vs. water) and some demands are subject to expiration while others may be carried over (e.g. food vs. clothing). The demands usually overwhelm the capacity of the distribution network. The demand information might not be complete and accurate at the beginning but it is expected to improve over time. Different aid organizations may employ their unique supply chain structure that governs the types of facilities to be used and the relationships among components of the chain. For example FEMA has its own supply chain structure for disaster response which is previously introduced in section 1.4. FEMA has distinguished 7 layers of facilities in its logistics chain. First 3 layers are permanent facilities to store and ship the relief items while the next 4 layers are temporary transfer facilities that their numbers and locations will be chosen during the response phase. 51 During the initial response time it is also necessary to set up temporary transfer facilities to receive, arrange, and ship the relief commodities through the distribution network. In risk mitigation studies for disasters, possible sites where these facilities can be situated are specified. Logistics coordination in disasters involves the selection of sites that result in the maximum coverage of affected areas and the minimum delays for supply delivery operations. Usually the number of these temporary facilities is limited because of the equipment and personnel constraints. Each facility in the chain is subject to some capacity constraints. Various capacities are defined for operations such as sending, receiving, and storing commodities. These capacities can be different for each facility and are determined based on the type, size and layout of that facility. Also the availability of personnel and equipment may influence the capacities. In general, the capacity constraints can be defined in terms of the weight or volume of the commodities as well as in terms of the numbers of the vehicles that are sent, received, or parked at the facility at a certain time period. These are two different aspects and it is recommended to consider both capacities for each facility. The transportation capacity is usually very limited in early hours or days after a disaster. It is very critical to assign the available fleet to the best possible use at any time. There is usually a shortage of vehicles in emergency operations so the model must keep track of the empty vehciles in order to assign them to new missions after each delivery. More than one transportation mode may be hired to facilitate emergency response logistics. Consequently, the coordination and 52 cooperation between transportation modes are necessary for managing the response operations and providing a seamless flow of relief commodities toward the aid recipients. The intermodal transfer of commodities is expected to happen in specific facilities but may be subject to some capacity constraints and transfer delays. Vehicle routing and scheduling during the disaster response is also very important. A large number of vehicles might be used in response to large scale disasters. The model should be able to keep track of routings for each individual vehicle. Also, it is required to have a detailed schedule for pickup and delivery of relief commodities by each vehicle in each transportation mode. Nonetheless, the vehicle routing in disaster situations are quite different from conventional vehicle routings. The vehicles do not need to form a tour and return to the assigned depot, but they might be assigned to a new path at any time. They are expected to perform mixed pickup and delivery of multiple items between different nodes of the network as the supplies and demands arise over time. The disaster area is a dynamic environment and emergency logistics are very time sensitive operations. The disaster might still be evolving when the response operations start. Also the lack of vital information about available infrastructure, supplies, and demands in the initial periods after the disaster may complicate this dynamic environment even more. The high stakes of life-or-death for disaster victims urge the needs for higher levels of accuracy and tractability. Despite all the necessary preparedness and planning at strategic level, dealing with the problem as operational level is very important. Modeling and 53 optimization at operation level is the only approach to capture the realities of the time sensitive emergency response operations. The other important issue is considering equity and fairness among aid recipients. Based on the geographical dispersion of victims and availability of resources over time and space, it is easy to favor the demands of one group of victims over another. Even though some variations are inevitable, the ideal pattern is to distribute the help items evenly and fairly among the victims. The models and procedures with general objective functions are prone to ignore the equity and level of service requirements in order to get a better numerical solution. It is very important to realize the need for procedures and constraints that prevent any sort of discrimination among victims, as much as possible. The equity constraint between populations can be defined over time, and over commodities. It is not appropriate to satisfy all the demands of one group in early stages while the other group of victims does not receive any help until very later times. It is more acceptable to fairly distribute the available relief items among all recipients even though it might not be enough for everyone at the current instance of time. The relief operations will continue over time as more resources are expected to become available. The equity over commodities is also important. For example, it is not acceptable to send all the available water to one group of victims and send all the available meals to another group. It is expected to fairly share the limited resources of transportation capacity and disaster relief commodities. 54 3.2 Time-Space Network A physical network is converted into a time-space network to account for the dynamic decision process. In the context of the problem in this research, nodes in the time-space network represent the physical locations of the supply, demand and transfer points for each mode and over time, while the arcs represent the connecting routes between these points. Each node in the physical network is represented by the number of mode types at each time period of the planning horizon. In a sense the time-space network in this context can be thought of as an overlay of several physical networks, one for each transportation mode, which are represented over time. These overlaid networks are connected to each other by the transfer links which make it possible for the commodities to be transferred between modes. There are three types of traffic flow on the physical network. The first type is the routing traffic that moves from one node to another node by a certain type of mode. The second type is the transfer traffic that changes mode type from one mode to another mode at a certain node. The third type is the supply or demand carry-over that is carried over to the next time period at a certain node. The duration of one time period should be based on the link travel time for each mode. It must be small enough so that the amount of slack time on the routing links is not excessive. However, the planning horizon should not be too short in order for the time-space network to be meaningful. Also, it should not be too long as it will increase the dimension of the time-space network and make the problem very difficult to solve. 55 The movements of commodities and personnel on a physical network over time are represented by the links in the time-space network. Routing Links represent the physical movement of commodities in space. Transfer Links represent the transfer of traffic between the available modes. Finally, Supply or Demand Carry-Over Links represent the commodity supply or demand carry over from one period to the next. Figure 3.1 shows a physical network that has 4 nodes, 5 two-way arcs, and 2 modes. Node A represents the origin and nodes C and D denote the destinations. The travel time over the arc in each mode type is shown in terms of time periods. Figure 3.2 shows the time-space network generated from Figure 3.1 with 6 time units in the planning horizon. The length of one time period is assumed to be one time unit. In Figure 3.2, transfer time is assumed to be one time period. The carry- over links that are created at node A and B represent the supply carry-over links. On the other hand, the carry-over links that are shown at node C and D denote the demand carry-over links. Figure 3.1 A sample physical network A C B D (1,2) (1,2) (2,4) (2,3) (1,2) Origin Destination Destination Transfer 56 Figure 3.2 A sample time-space network 3.3. Modeling Approach A mathematical framework is suggested to model the supply chain operations during emergency response similar to the problem description in section 3.1. The main characteristics of the modeling approach can be summarized as follow: ? Operational Level: to capture time sensitive details of the emergency response operations, the problem is formulated at operational level. ? FEMA Structure: the proposed model is in compliance with FEMA?s 7-layer supply chain structure. Mode M Time period 0 1 2 3 4 5 6 Routing Flow Transfer Flow Carry-over Flow A B C D A B C D Mode N 57 ? Time-Space Network: to account for the dynamic decision process, the physical network must be converted to a time-space network. The nodes of this network represent the facilities in the FEMA?s structure. The links consist of existing physical links, delay or storage links, and intermodal transfer links. ? Facility Location: the optimal locations to establish temporary facilities are selected from a set of potential sites. The maximum number of each facility type and their locations are dynamic and can change over time as the relief operations proceed. ? Facility Capacity: each facility has maximum capacities for sending, receiving, and storing commodities as well as vehicles. ? Demand: the demand is multi-commodity and usually overwhelms the capacity of the distribution network. Specific decision variables are defined that keep track of unsatisfied demand at each demand point for each commodity and during all time periods. ? Supply: similar to the demand, the supply is multi-commodity and may come from various sources. The problem is formulated as a variation of multi- commodity network flow problem. ? Multi-modal: since more than one mode of transportation may be hired in the emergency response logistics, the problem is a variation of multi-modal network flow problem. ? Vehicle Routing: in order to model the complicated routing and delivery operations in disaster response, the vehicles are treated as flow of integer 58 commodities over the time-space network. This results in a mixed integer multi-commodity formulation. ? Network Capacity: a set of constraints is used to link the relief commodities with the vehicles. As a result, the flow of commodities is only possible when accompanied by enough vehicle capacity for that specific link and time. ? Integrated Model: all decisions of facility location, supply delivery, and vehicle routing, are interrelated. Our approach provides an integrated model to find the global solution for this problem. ? Equity: equity and fairness among disaster victims is modeled through a set of constraints that enforce a minimum level-of-service for each victim. The equity is required for each relief item and over all time periods. ? Objective Function: the objective of this model is to minimize the pain and suffering of the disaster victims. It is formulated as minimizing the total unsatisfied demand summarized for all victims, for all relief items, and during all time periods. 3.4 Assumptions 1- It is assumed that the following information is available and given: ? Demands: commodity types, demand locations, demand amounts ? Supply: commodity types, supply locations, supply amounts ? Permanent Facilities: types, locations, capacities ? Temporary Facilities: set of potential sites for each type, capacities of each type 59 ? Network: link-node incidence matrix for each transportation mode ? Vehicles: number of vehicles available for each mode and their initial locations, capacity of each vehicle ? Travel Times: travel time on each link for each transportation mode. 2- Because the model is at the operational level, it is assumed that the problem is deterministic. The required information is estimated or known at the beginning of the operations by local or federal authorities. The model can adapt to the new information as the circumstances evolves over time and real-time information becomes available. 3- Supply Chain Structure: ? It is assumed that the flow of commodities between each two nodes is possible only if it is in compliance with FEMA?s structure shown in Figure 1.5. For example, the supply from LC cannot be sent directly to SSA. It should be sent to MOBs or FOSAs first. ? It is assumed that for the empty vehicles, a direct link exists that connects each pair of nodes. For example, if a vehicle delivers all of its supply at a POD, it can directly go to any other node of the network to pick up new supplies. 4- Finding the number and locations of the points of distribution (PODs) is not considered in this study. It is assumed that PODs are established by local authorities. As a result, the location and amount of demands at each POD is a given data in this model. 60 3.5 Mathematical model In this section initially the notations and required parameters for the formulation are introduced. After that, the decision variables of the mathematical model are defined. Then the objective function formulation is presented followed by introduction and formulation of the problem constraints. 3.5.1 Notations N = Set of all nodes. Nji ?, are indices LC = Set of logistic center sites CSS = Set of commercial storage sites VEN = Set of commodity vendor sites MOB = Set of potential sites for mobilization centers FOSA = Set of potential sites for federal operational staging areas SSA = Set of potential sites for state staging areas POD = Set of points of distribution (demand nodes) U = Set of supply nodes and transshipment nodes (LC, VEN, CSS, MOB, FOSA, SSA) V = Set of permanent facilities (LC, CSS, VEN) W = Set of potential sites for all temporary facilities (MOB, FOSA, SSA) C = Set of commodities, Cc ? is an index M = Set of transportation modes, Mm ? is an index T = Time horizon of response operations. Ttt ??, are indices 61 3.5.2 Parameters Supply and Demand c itSup = Amount of exogenous supply of commodity type c in node i at time t c itDem = Amount of exogenous demand of commodity type c in node i at time t m itAV = Number of vehicles of mode m added to the network in node i at time t, negative if vehicles removed c itRU = Relative urgency of one unit of commodity c, in node i at time t Number of Facilities tMOB max = Maximum number of mobilization centers at time t tFOSA max = Maximum number of federal operational staging areas at time t tSSA max = Maximum number of state staging areas at time t Facility Capacity m itUcap = Unloading capacity for the facility in node i for mode m at time t itScap = Storage capacity for the facility in node i at time t m itLcap = Loading capacity for the facility in node i for mode m at time t m itVRcap = Maximum number of mode m vehicles that can be received at the facility in node i at time t m itVPcap = Maximum number of mode m vehicles that can be parked (carried over) at the facility in node i from time t to time t + 1 m itVScap = Maximum number of mode m vehicles that can be sent out from the facility in node i at time t 62 Vehicle Capacity mcap = Loading capacity of vehicles of mode m cw = Unit weight of commodity c Transportation ijmt = Travel time from node i to node j for vehicles of mode m mmK ? = Time required to transfer commodities from mode m to mode m? 3.5.3 Decision Variables Location Problem t iLoc = 1 if temporary facility of appropriate type is located at potential site i, at time t; equal to 0 otherwise. The temporary facility will be a mobilization center if MOBi ? , a federal operational staging area if FOSAi ? , and a state staging area if SSAi ? . Commodity and Vehicle Flow cm ijtX = Flow of commodity type c shipped from node i to node j by mode m at time t m ijtY = Flow of vehicles of mode m from node i to node j at time t c itCX = Amount of commodity type c in node i which is carried over from time period t to t + 1 m itCY = Number of vehicles of mode m in node i which is carried over from time period t to t + 1 63 mcm itXT ? = Amount of commodity type c in node i which is transferred from mode m to mode m? at time t c itUD = Amount of unsatisfied demand of commodity type c in node i at time t 3.5.4 Objective Function Minimize ??? ? ? Vi t c c it c it UDRU (3.1) The objective function in equation (3.1) minimizes the total amount of weighted unsatisfied demand over all commodities, times, and demand points. c itRU is the relative urgency associated with each commodity, time, and demand point. If there is any desire to consider a commodity being more important than others at any time or for any demand point, citRU can enforce that desire. Higher value of citRU translates into higher urgencies. If all the commodities happen to be of the same importance, citRU can be set equal to 1. 3.5.5 Constraints Commodity Flow Constraints Supply nodes and Transfer nodes: tmcUiCXXTX SupCXXTX c it m mcm it j cm ijt c it c ti m mmc kti j cm ttji mmjim ,,, )1()()( ??++= +++ ?? ?? ? ? ? ? ? ?? ? (3.2) Demand nodes: tcPODiUDDemUDX c ticitcit m j cm ttji jim ,,)1()( ??+=+ ???? (3.3) 64 Equations (3.2) and (3.3) enforce the conservation of the flow for all commodities and modes at all nodes and time periods. Equation (3.2) requires that for supply nodes and transfer nodes, the sum of the flows entering each node plus exogenous supply should be equal to the sum of the flows that leave the same node. Equation (3.3) shows that the total flow entering each demand node plus the unsatisfied demand is equal to the exogenous demand at that node plus the unsatisfied demand from the previous time period. Vehicular Flow Constraints tmNiCYYAVCYY mit j m ijt m it m ti j m ttji jim ,,)1()( ??+=++ ?? ?? (3.4) Equation (3.4) represents the conservation of flow for the vehicles. At any node i and time period t, total number of available vehicles of mode m is equal to the number of vehicles of mode m that left node j for node i at time ijmtt ? , plus the number of vehicles that were carried over from the previous time period, plus the number of vehicles that are added or removed to the fleet at that time. These vehicles are either sent out of the node or carried over to the next time period. Linkage between Commodities and Vehicles tmNiXwYCap c cm ijtc m ijtm ,,???? ? (3.5) Constraint (3.5) makes sure that commodities are not sent out of a node unless a number of vehicles with enough capacity are available at that node. Facility Capacities for Permanent Facilities tmViLcapX mit c j cm ijt ,,????? (3.6) 65 tmLCiUcapX mit c j cm ttji jim ,,)( ????? ? (3.7) tViScapSupCXX it c c it c c ti m c j cm ttji jim ,)1()( ???++ ????? ?? (3.8) tmViVScapY mit j m ijt ,,???? (3.9) tmViVRcapAVY mitmit j m ttji jim ,,)( ???+? ? (3.10) tmViVPcapCYAVY mitmtimit j m ttji jim ,,)1()( ???++ ??? (3.11) Equations (3.6), (3.7), and (3.8) are the maximum capacity for loading, unloading, and storage of commodities at permanent facilities. Equations (3.9), (3.10), and (3.11) require the maximum number of vehicles that are sent, received, and parked at each facility to be less than the relevant capacities. Facility Location and Capacities for Temporary Facilities tmWiLocLcapX timit c j cm ijt ,,?????? (3.12) tmWiLocUcapX timit c j cm ttji jim ,,)( ?????? ? (3.13) tWiLocScapSup CXX t iit c c it c c ti m c j cm ttji jim , )1()( ????+ + ? ???? ?? (3.14) tmWiLocVRcapAVY timitmit j m ttji jim ,,)( ????+? ? (3.15) tmWiLocVPcapCYAVY timitmtimit j m ttji jim ,,)1()( ????++ ??? (3.16) tmWiLocVScapY timit j m ijt ,,????? (3.17) 66 tMOBiMOBLoc t i t i ,max ???? (3.18) tFOSAiFOSALoc t i t i ,max ???? (3.19) tSSAiSSALoc t i t i ,max ???? (3.20) Equations (3.12) through (3.14) enforce the loading, unloading, and storage capacity for the temporary facilities. If the facility is selected to be set up at potential site i, the respected capacity constraint is enforced. If it is decided not to set up the temporary facility at location i, the same constraints require that all the flows in and out of that node to be equal to zero. Equations (3.15) through (3.17) require the maximum number of vehicles that are sent, received, and parked at each temporary facility to be less than the relevant capacities. The numbers are zero if the facility is not selected for that node. Equations (3.18) through (3.20) oblige the maximum number of each temporary facility type to be limited by the maximum allowable numbers for that facility type during the chosen time periods. Capacities for PODs: tmPODiUcapX mit c j cm ttji jim ,,)( ????? ? (3.21) tmPODiVRcapY mit j m ttji jim ,,)( ???? ? (3.22) tmPODiVPcapCYY mitmti j m ttji jim ,,)1()( ???+ ??? (3.23) 67 Equation (3.21) enforces the commodity unloading capacity at points of distribution. Equation (3.22) and (3.23) represent the vehicle receiving and vehicle parking capacities for each point of distribution. Equity Constraint: tcPODiDem X t c ti t m j cm ttji jim ,,min )( ???? ??? ? ? ? ?? a (3.24) tcPODiX X c t m j cm ttji t m j cm ttji jim jim ,,min )( )( ??????? ??? ? ?? ? ?? b (3.25) tPODiX X i c t m j cm ttji c t m j cm ttji jim jim ,min )( )( ???????? ???? ? ?? ? ?? g (3.26) Equation (3.24) enforces a minimum percentage of total demand for a specific commodity c, to be satisfied by the time period t. It might not be always possible to deliver the required amount by time t; in that case, this constraint makes the optimization problem infeasible. Equation (3.25) requires that from all commodities being delivered to node i by time t, at least minb percent to be commodity c. Equation (3.26) ensures that sum of total commodities delivered at point i to be more than a minimum percentage of all the commodities that are being delivered among all demand points. Non-negativity and Integrality: cm ijtX , c itCX , mcm itXT ? , 0?c itUD Real-valued variables 68 m ijtY , 0? m itCY General integer variables )1,0(?tiLOC Binary integer variables 3.6. Summary The proposed mathematical model in this chapter can be summarized as follows: Minimize Total weighted unsatisfied demand Subject to: Commodity flow constraints Vehicular flow constraints Constraints that link commodities and vehicles Facilities location constraints Facility capacities constraints Equity (recipients/commodities) constraints Non-negativity and integrality constraints 69 Chapter 4: Preliminary Numerical Study In this chapter, a set of preliminary numerical experiments are conducted to evaluate the features of the proposed formulation. At this stage, it is tried to keep the problem size manageable so it can be solvable by commercial solver and the results can be analyzed easier. Nevertheless, this problem instance still fully represents all the elements of the proposed model. This experimental study is compliant with FEMA?s special supply chain structure. The distances, locations, supplies, demands, capacities and other aspects of this numerical experiment are designed to be comparable to the real-world-size problems. 4.1 Design of the sample problems The numerical problem in this chapter is an imaginary scenario where a natural disaster such as a hurricane strikes the southern coast of the United States. It is assumed that two separate regions, one in Mississippi and one in Louisiana, are affected. The disaster area in Mississippi is spread along the coast while the disaster area in Louisiana is more inland and has a rectangular shape. Figures 4.1 and 4.2 show the affected disaster areas. 70 Figure 4.1 General Disaster Area of the Numerical Study Figure 4.2 Disaster Areas in two affected States Mississippi Louisiana 71 4.1.1 NETWORK For the numerical study, it is assumed that only the Atlanta logistics center (LC) is used. One commercial storage site (CSS) in Charlotte, North Carolina and one vendor (VEN) in Nashville, Tennessee are also used to store the relief items. For temporary facilities at federal level, four potential sites for mobilization centers (MOB) are suggested. There are also four potential sites for federal operational staging areas (FOSA). These facilities are able to send supplies to both disaster areas. At the state level, a total of 10 potential sites for state staging areas (SSA) are suggested. Four potential SSA are planned to serve the disaster area in Mississippi and six potential SSA are suggested for Louisiana. The initial post-disaster surveys estimate that approximately 20,000 people are affected and twenty points of distribution (POD) are needed to serve this population. Eight PODs are selected for Mississippi area and twelve PODs will serve the victims in Louisiana. Table 4.1 summarizes the list of facilities in the distribution network. For this numerical experiment, there are a total of 41 permanent and temporary facilities in the network. Figures 4.3, 4.4, and 4.5 illustrate the locations of these facilities on the map. 72 Table 4.1 List of Facilities in the Distribution Network Node Facility TYPE Location Latitude Longitude 1 LC Atlanta, GA 33?44'6.59"N 84?23'45.13"W 2 CSS Charlotte, NC 35?13'47.00"N 80?50'36.54"W 3 VEN Nashville, TN 36?11'9.34"N 86?43'25.24"W 4 MOB Montgomery, AL 32?22'3.90"N 86?18'6.88"W 5 MOB Jackson, MS 32?18'20.21"N 90?10'7.65"W 6 MOB Shreveport, LA 32?30'44.00"N 93?44'25.76"W 7 MOB Beaumont, TX 30? 4'47.76"N 94? 6'2.57"W 8 FOSA Mobile, AL 30?41'20.63"N 88? 2'44.56"W 9 FOSA Hattiesburg, MS 31?18'16.67"N 89?18'41.34"W 10 FOSA Baton Rouge, LA 30?26'49.07"N 91?11'4.33"W 11 FOSA Lafayette, LA 30?12'39.24"N 92? 0'36.65"W 12 SSA Moss Point, MS 30?25'36.88"N 88?31'20.06"W 13 SSA Gulf Hills, MS 30?26'14.86"N 88?48'52.52"W 14 SSA Wool Market, MS 30?28'4.60"N 88?59'49.49"W 15 SSA Diamond Head, MS 30?22'48.38"N 89?22'32.34"W 16 SSA Boutte, LA 29?54'5.23"N 90?23'28.72"W 17 SSA South Vacherie, LA 29?54'40.81"N 90?43'44.11"W 18 SSA Supreme, LA 29?52'2.73"N 90?59'4.48"W 19 SSA Pierre Part, LA 29?57'19.71"N 91?12'45.39"W 20 SSA Berwick, LA 29?42'3.16"N 91?13'51.50"W 21 SSA Franklin, LA 29?47'17.49"N 91?30'33.94"W 22 POD Pascagoula, MS 30?21'54.42"N 88?32'54.99"W 23 POD Gautier, MS 30?23'26.03"N 88?38'44.36"W 24 POD Gulf Park, MS 30?22'45.27"N 88?45'32.84"W 25 POD Ocean Springs, MS 30?24'39.92"N 88?47'7.53"W 26 POD Biloxi, MS 30?24'27.58"N 88?55'59.03"W 27 POD Gulf Port, MS 30?21'57.06"N 89? 5'30.75"W 28 POD Long Beach, MS 30?20'24.34"N 89?11'1.03"W 29 POD Pass Christian, MS 30?19'33.94"N 89?14'57.81"W 30 POD Lock Port, LA 29?38'22.61"N 90?32'14.66"W 31 POD Mathews, LA 29?41'38.04"N 90?33'6.94"W 32 POD Raceland, LA 29?43'19.20"N 90?35'17.82"W 33 POD Houma, LA 29?35'13.92"N 90?42'15.67"W 34 POD Bayou Cane, LA 29?37'29.72"N 90?45'3.30"W 35 POD Gray, LA 29?40'45.88"N 90?47'0.88"W 36 POD Shriever, LA 29?44'25.98"N 90?49'50.30"W 37 POD Tibodaux, LA 29?47'48.50"N 90?49'7.77"W 38 POD Amelia, LA 29?40'16.24"N 91? 6'15.78"W 39 POD Morgan City, LA 29?42'9.13"N 91?11'25.60"W 40 POD Bayou Vista, LA 29?41'28.15"N 91?16'13.42"W 41 POD Patterson, LA 29?41'23.98"N 91?18'33.41"W 73 Figure 4.3 Map of the federal level facilities Figure 4.4 Map of the state level facilities in Mississippi LC(1) CSS(2) VEN(3) MOB(4)MOB(5)MOB(6) MOB(7) FOSA(8) FOSA(9) FOSA(10)FOSA(11) SSA(12) SSA(13) SSA(14) SSA(15) (22) (23)(24) (25)(26) (27) (28)(29) 74 Figure 4.5 Map of the state level facilities in Louisiana 4.1.2 Supply and Demand There are several commodities that need to be distributed among the disaster victims. The type and amount of each commodity depends on many factors such as type of disaster, level of destruction, weather conditions, etc. Table 4.2 suggests a list of required items and the amount per day per survivor. Adding up the last column of Table 4.2, it can be seen that for each survivor a total of about 30 ft3 of relief items per day are required. For the sake of simplicity, it is assumed that only 2 types of commodities (commodity 1 and commodity 2) are required in this numerical experiment. However, to preserve the scale of demands, the total amount per each survivor is kept at 30 ft3 per day. It is also assumed that survivors in disaster zone 1 (Mississippi), need 20 ft3 of commodity 1 and 10 ft3 of commodity 2, per day. On the other hand, survivors in disaster zone 2 (Louisiana), assumed to need 10 ft3 of SSA(16) SSA(17) SSA(18) SSA(19) SSA(21) SSA(20) (30) (31) (32) (33) (34) (35) (36) (37) (38) (39)(40)(41) 75 commodity 1 and 20 ft3 of commodity 2, per day. This will provide the opportunity to analyze the effects of different demand types in the results of the model. Table 4.2 List of Required Items for Survivors of a Disaster Item Quantity Survivors served dimensions (ft) Volume (ft3) Requirement per survivor (ft3) L W H Water (drinking) 1 gallon 1 1 1 1 1 3 Water (non-potable) 1 gallon 1 1 1 1 1 3 Meals (MREs) 3 meals 1 1 1 1.5 1.5 4.5 Portable shelter 1 shelter 4 6 2 1.5 18 4.5 Basic medical kit 1 kit 3 1 1 1 1 0.333 Cot 1 cot 2 3 2 1 6 3 Blanket 1 blanket 1 2 2 0.5 2 2 Tarp 1 tarp 3 3 3 1 9 3 Ice 1 gallon 10 1 1 1 1 0.1 Baby supplies 1 box 5 1 1 1 1 0.2 Generator 1 generator 500 8 8 6 384 0.768 Clothing 1 bag 1 2 2 1 4 4 Plywood 2 sheets 3 4 8 0.1 3.2 2.133 Nails 1 box 3 1 1 1 1 0.333 Supply sources are the logistics center, the commercial storage site, and the vendor. It is assumed that 40% of total supply is stored at LC, 20% at CSS, and 40% at the vendor site. Total demand for 20,000 survivors will be 600,000 ft3 per day. The demand for commodity 1 is 280,000 ft3 per day and the demand for commodity 2 is 320,000 ft3 per day. For this problem, it is assumed that supplies for one day are available and are stored at the three supply sources prior to the start of the operation. 76 4.1.3 Transportation For this problem, only one transportation mode is used which is trucking. The common vehicle is a 53ft trailer truck which has the volume capacity of approximately 6000 ft3. For the base case, 100 trucks are available at the beginning of the operations. Initially, 40 trucks are located at LC and 30 trucks are present at CSS and VEN sites, each. 4.1.4 Network links and travel times There are 2 types of flows in this problem, flow of commodities and flow of vehicles. The commodity flows must comply with the hierarchical structure of FEMA explained in section 1.4. For example, supplies from a VEN can only be sent to LC, or supply from LC can be sent to all MOBs and FOSAs. Supplies in MOBs can be sent to other MOBs or to FOSAs. Supplies from FOSAs can be sent to other FOSAs and to SSAs, as long as it remains in the same state. Supplies received at each SSA can be sent to other SSAs in the same State or must be delivered to PODs of that State. The flow of vehicles in the network is much less restricted compared to commodity flows. It is assumed that there is a link between each pair of nodes in the network. Basically, empty vehicles are free to travel between each two nodes of the network without the need to visit any intermediate nodes. As a result, when a vehicle is carrying supplies, it must follow the more restricted hierarchical network of FEMA. But when the vehicle unloads all its supply, either at intermediate nodes or final PODs, it is free to go to any other node in the network to pick up supplies and start a new round of deliveries. 77 Link travel time functions for the proposed formulation can be completely arbitrary. The formulation is capable of dealing with time-variable travel times as well as fixed travel times. For this numerical study, the travel distance between any two nodes of the network is assumed to be equal to their Euclidian distance. The travel speed is assumed to be fixed for all the vehicles on the federal level network (between LC, CSS, VEN, MOBs, and FOSA) and to be equal to 50 miles per hour. However, for the state level networks (between FOSAs, SSAs, and PODs) the travel speed is assumed to be 40 miles per hour. 4.1.5 Time Scale Selection of appropriate time step is a very important factor that can affect the performance of time-space networks dramatically. For each time period in the planning horizon, one layer of physical network will be added to the problem. This makes the problem size grow extremely fast with the number of time steps in the planning horizon. For example if the planning horizon is only 1 day, with the choice of time step t = 5 minutes, it will be 24 * 60 / 5 = 288 layers of the network. So to keep the problem at a reasonable size, it is favorable to have longer time steps. On the other hand, shorter time steps will improve the accuracy of modeling the emergency response operations. For example if the time step is 1 hour, it is possible to model the state of the system only at every hour and not at the times in between. So from the accuracy perspective, it is favorable to have shorter time steps. 78 The other important issue in determining the time-step in this problem is the issue of dealing with very long and very short links. At the federal level network, nodes are usually far from each other and the links can range from some hundred miles to a few thousand miles. The travel time on those links with ground transportation can range from a few hours to up to one day or more. However, the nodes at the lower levels in the State networks can be very close to each other. It is very common to have PODs that are only a few miles apart. In this case, link travel times can be in the order of minutes. Figure 4.6 better shows the issue of scale in this problem in disaster area map. Figure 4.6 The issue of the scale in the disaster area It is a difficult challenge to select a time-step that is suitable for very short links and very long links, at the same time. A very short time-step is necessary to model the short links even though it will increase the problem size very quickly. But the main issue is the sensitivity of travel times to the selected time-step. If a very short time-step is chosen, say 1 minute, it might be good for short links but 79 the travel times on very long links will not be sensitive to that. It is very difficult, if not impossible, to predict the travel time between two nodes that are a thousand miles apart, with accuracy of 1 minute. For those links the 1-hour unit or 30- minute unit is more meaningful. 4.1.6 Geographical Decomposition To deal with the issue of scale, a geographical decomposition method is proposed. The nodes at federal level (LC, CSS, VEN, MOB, FOSA) will be in one subset and the nodes at each state (FOSA, SSA, POD) will form another subset. Since the travel times between nodes in federal level network are usually long, it is possible to use a large time-step for them. Using similar argument, the State level nodes and links can be modeled with a shorter time-step. Figure 4.7 shows this decomposition. Figure 4.7 Geographical decomposition to implement two time-steps VEN LC CSS MOB MOB FOSA FOSA SSA SSA SSA SSA POD POD POD POD POD POD fed St1 St2 80 Now the important issue is how to connect these separate time-space networks. Luckily, the special structure of FEMA?s supply chain offers the candidates. Federal operational staging areas (FOSA) are the one and only interface between flow of commodities among the federal level facilities and the designated state level facilities. We take advantage of this opportunity and select the FOSAs as transfer terminals between the networks with different time steps. For this numerical study, time-step for federal zone, 1t , is chosen to be 30 minutes and time-step for state level zones, 2t , is selected to be 5 minutes. The travel times for this study are calculated based on the distance and a fixed average travel speed explained earlier. So based on the newly defined time steps of 1t and 2t , travel times of federal zone links are being rounded to the nearest 30 minute interval and the travel times of state level zone links are being rounded to the nearest 5 minute. The way FOSA nodes connect two sub-networks with different time steps is shown in Figure 4.8. This graph indicates that the arcs entering FOSA from federal network or leaving the FOSA toward the federal network can exist only at 1t =30-minute intervals. But the arcs that connect FOSA to state level facilities exist for every 2t =5-minute interval. The implication is that the downward flows (from the federal network to the state network) entering a given FOSA can leave that FOSA at any 5-minute period after that. However, the upward flows (from the state networks to the federal network) that enter a FOSA at any time other than 30-minute intervals need to wait at the FOSA until the first available 30- minute interval. 81 Figure 4.8 Connecting two time-Space networks with different time 4.2 Generating formulation for commercial solver The numerical experiment introduced previously in section 4.1 is a fairly large mixed integer program with real valued as well as general integer and binary variables. Because of the large number of variables and constraints in this problem, computer programming is required to handle the input and output data. A customized program is coded in the Microsoft Visual Studio environment to generate the mathematical formulation for each problem instance. The program read the input data from the prepared data files as well as the coded user interface to generate each problem instance. Then the mathematical formulation for each problem instance is generated and written to a text output file. FOSA State-Level Nodes Federal-Level Nodes Time period 0 1 2 3 4 5 6 7 8 9 10 11 12 t2 t1 Downward Flow Upward Flow Carry-over Flow 82 At this stage, the numerical sample problems are solved with ILOG CPLEX (2006) commercial solver. CPLEX is a commercial optimization package from ILOG company that solves mathematical formulations in the forms of linear programs (LP), integer programs (IP), and quadratic programs (QP). CPLEX reads the generated formulations from text files, after optimization the results are written to text files as well. Another customized program is coded to extract the results from the output file and generate the required performance measures and generate charts and graphs. 4.3 Numerical results and analysis To better evaluate the characteristics of the proposed model, 10 numerical case studies are generated. All the case studies are based on the described imaginary scenario with variations in the subset of enforced constraints and some parameter values. Table 4.3 describes the considered case studies. In general, the case studies in tables 4.3 start from simple and become more complicated toward the end. For example, the first case study only considers the conservation of flow and vehicle capacity constraints. Other constraints are gradually added to the formulation in the other case studies up to case 7 which has the largest number of constraint types for a one day operation. First 7 case studies consider only 1 day of operations while in the last 3 cases 2 days of operations are formulated. Table 4.4 summarizes the optimization results for all 10 case studies. case- 1 is the ?base case? with only conservation of flow constraint and vehicle capacity constraints modeled for 1 day of operations. The solver found the optimal solution in approximately 4 minutes. Figure 4.9 shows the percent of unsatisfied demand 83 for all victims over time. The first delivery to the nearest demand point took about 7 hours. Fifty percent of the total demand was satisfied after 11 hrs and 40 minutes. The last demand was served after 21 hours and 40 minutes. Table 4.3 Descriptions of the numerical Case Studies Case Constraints Used Details No. of Variables Constraints File Size (Kb) Real Integer 1 Flow Conservation + Vehicle Capacity 1 day 100 Trucks 133275 157972 81,891 13,331 2 Flow Conservation + Vehicle Capacity 1 day 200 trucks 133275 157972 81,891 13,331 3 Flow + Vehicle Capacity + Facility Capacity 1 day 100 Trucks 133275 157972 87,094 15,846 4 Flow + Facility Location (2,2,5)* + Facility Cap 1 day 100 Trucks 133275 157972 87,094 15,846 5 Flow + Facility Location (2,2,2) + Facility Cap 1 day 100 Trucks 133275 157972 87,094 15,846 6 Flow + Facility Capacity Const.+ Equity-1 Const 1 day 100 Trucks 133275 157972 87,174 17,214 7 Flow + Facility Location & Capacity + Equity-1,2,3 1 day 100 Trucks 133275 157972 87,294 61,084 8 Flow + Vehicle Cap, day by day Supply 2 days 100 Trucks 265995 315316 163443 27,439 9 Flow + Facility Location & Capacity, day by day 2 days 100 Trucks 265995 315316 173,878 32,673 10 Flow + Capacity + location (2,2,5) , 2 day supply available 2 days 100 Trucks 265995 315316 173,878 32,673 * Facility location with maximum number of (MOB, FOSA, SSA) 84 Table 4.4 Summary of the optimization results for the preliminary experiment Case Number Objective Value (E+7) Last UD (hr:min) Temp. Facilities Root Sol. Time (s) Iterations CPU Time (sec)? 1 9.0798 21:40 (4,4,10) 33.89 14,957 230 2 8.6118 15:10 (4,4,10) 10.36 5,502 20 3 10.412 22:05 (4,4,10) 42.73 18,642 778 4 10.412 22:05 (2,2,5) 33.59 17,308 945 5 10.978 24:00? (2,2,2) 204.19 205,588 5575 6 10.439 21:50 (4,4,10) 42.22 5,810,980 45856* 7 10.417 22:05 (2,2,5) 63.09 7,888,315 81642* 8 17.985 39:10 (4,4,10) 786.34 63,960 4779 9 20.859 44:45 (2,2,5) 2450.91 408,351 14635 10 18.921 48:00? (2,2,5) 10117.11 2,963,071 231035 * The solver stopped prematurely with ?out of memory? error message. ? The relief operations were not finished by the end of planning horizon. ? On a 3.0 GHz Intel Pentium CPU with 2.0 GB RAM Case-2 is similar to case-1 but the only difference is that there are 200 trucks available in case-2 versus 100 trucks in case-1. Even though the number of vehicles was increased, the optimal solution was found in only 20 seconds. As it can be seen in Table 4.3, the size of the formulation (number of variables and constraints) for case-2 is equal to case-1 and this is one of the important advantages of current formulation. Since this formulation treats the vehicles as commodities, the number of available vehicles appears only as a right-hand-side parameter and does not have an effect on the problem size. Figure 4.10 shows the 85 percent of unsatisfied demand over time for case-2 at optimality. Since there were enough vehicles available at the beginning, the vehicles did not need to return to the sources to pickup supplies once they had left. As a result, the delivery operations were completed after only 15 hours and 10 minutes. Figure 4.9 Percent of unsatisfied demand over time for case 1 Figure 4.10 Percent of unsatisfied demand over time for case 2 0 10 20 30 40 50 60 70 80 90 100 0 2 4 6 8 10 12 14 16 18 20 22 24 Time (hr) Un sa tis fie d D em and (%) 0 10 20 30 40 50 60 70 80 90 100 0 2 4 6 8 10 12 14 16 18 20 22 24 Time (hr) Un sa tis fie d D em and (%) 86 Case 3 is similar to the base-case with addition of loading, unloading, and storage capacities for all facilities. In this case, there is no limitation on the maximum number of temporary facilities and all the potential sites can be active. Figure 4.11 shows the variation of unsatisfied demand for case-3. The addition of facility capacities prevented the shipment and delivery of large quantities of supplies. Instead, the relief commodities are delivered more uniformly over time compared to case-1 and Figure 4.9. Consequently, the objective function value was higher and the operations lasted for 22 hr and 5 minutes, 25 minutes more than case-1. The running time was also increased to about 13 minutes to find the optimal solution. In case 4 we limited the maximum number of temporary facilities (MOB, FOSA, SSA) to (2, 2, 5) plus the constraints of case-3. It took the solver about 16 minutes to find the optimal solution which is 3 minutes more than case-3. However, the objective function value at optimality was the same for case-3 and case-4. This implies that although we limited the number of temporary facilities to (2,2,5); it was still possible to run the operations through limited number of facilities and achieve the same final results. Comparing Figure 4.12 with 4.11 shows that there were minor changes in the flow of commodities, but the final results are very similar. 87 Figure 4.11 Percent of unsatisfied demand over time for case 3 Figure 4.12 Percent of unsatisfied demand over time for case 4 In order to see the effect of even more limited numbers of facilities, we created case-5 with maximum number of temporary facilities as (2,2,2). Table 4.4 shows that the problem became much harder to solve. The running time jumped to 1.5 hours. The objective function was increased and more importantly, the 0 10 20 30 40 50 60 70 80 90 100 0 2 4 6 8 10 12 14 16 18 20 22 24 Time (hr) Un sa tis fie d D em and (%) 0 10 20 30 40 50 60 70 80 90 100 0 2 4 6 8 10 12 14 16 18 20 22 24 Time (hr) Un sa tis fie d D em an d ( %) 88 delivery operations could not be finished in 24 hours. Figure 4.13 shows that at the end of the 24 hours, there is still unsatisfied demand which is about 6% of the total demand. It indicates that unlike case-4, limiting the number of temporary facilities affected the operations and resulted in more delays and more unsatisfied demand. In case-6 we added the equity constraints to the problem for the first time. At this stage, only the 1st equity constrains (Equation 3.24) were considered in addition to conservation of flow and vehicle capacity constraints. Table 4.4 shows very interesting results. First of all, adding the equity-1 constraint made the problem much harder. After 13 hours of execution time and more than 5.8 million iterations, the solver still could not find the optimal solution. However, the best integer solution found is very close to the best MIP bound (25500 unsatisfied demands, 0.02% gap). Figure 4.13 Percent of unsatisfied demand over time for case 5 0 10 20 30 40 50 60 70 80 90 100 0 2 4 6 8 10 12 14 16 18 20 22 24 Time (hr) Un sa tis fie d D em and (%) 89 Finally in Case-7, all the constraints are considered. The constraints include conservation of flow for the commodities and vehicles, the linkage between commodities and vehicles and capacity of each vehicle, facility location with maximums of (2,2,5); loading, unloading and storage capacities for all facilities, and finally the 3 equity constraints (Equation 3.24, 3.25, 3.26). The full problem becomes very large and difficult to solve. After around 23 hours of CPU time and more than 7.8 million iterations, CPLEX solver stopped and it could not find the optimal solution. By the way, the best integer solution found is very close to the best MIP bound (30400 unsatisfied demands, 0.03% gap). Figure 4.14 shows the unsatisfied demand for the best integer solution found by the solver. Figure 4.14 Percent of unsatisfied demand over time for case 7 Another idea was to extend the relief operations duration from 1 day to 2 days and analyze its effect on the problem size and behavior. Case-8 through case-10 was created to test this idea. Again case-8 is the base case with conservation of flow and vehicle capacity constraints. Similar to other cases, it 0 10 20 30 40 50 60 70 80 90 100 0 2 4 6 8 10 12 14 16 18 20 22 24 Time (hr) Un sa tis fie d D em and (%) 90 was assumed that 100 trucks are available. The demands for the second day appear at the beginning of the second day and locations and quantities are similar to the demands of the first day. In case-8 it is assumed that supply for the second day arrives at the beginning of the second day of operations to the same source nodes as in day one. Table 4.4 shows that the solution time was around 80 minutes. Comparing 80 minutes for case-8 with 4 minutes of CPU time for case-1 shows the growth rate of problem size and difficulty with extending time horizon. In this case, the duration of operation is only doubled however the solution time is rapidly increased by a factor of 20. Figure 4.15 shows the variations in unsatisfied demand over time. The 1st day?s operations were finished in approximately 18 hours. As a result of the optimal distribution of empty trucks for the second day, the relief operations in the second day were over in only 15 hours and 10 minutes. There were no additional supplies available before the second day, but modeling the operations for 2-day provided the opportunity to be prepared and do a better job in the second day. In case-9, the facility location constraints with maximum of (2,2,5) and the loading, unloading, and storage capacity constraints were considered for 2 days of operations. Similar to the previous case, the supplies become available day by day. Table 4.4 shows that adding the capacity constraints has increased the objective function value for about 16% compared to case-8. The running time is also increased to more than 4 hours of CPU time. 91 Figure 4.15 Percent of unsatisfied demand over time for case 8 Figure 4.16 shows the results of the optimal solution for case-9. Because of the capacity constraints the flow of extra large amounts of commodities were prohibited. As a result, the demands are satisfied gradually over time and for both days, the operations took longer compared to case-8. It took 44 hours and 45 minutes to deliver all the supplies in case-9 compared to 39 hours and 10 minutes in case-8. The last case study in this numerical experiment is case-10. case-10 is similar to case-9 with the only variation that all the supplies for 2-days are assumed to be available at the beginning of the operations. The demands still appear at the start of each day and supplies cannot be delivered beforehand. The objective function of optimal solution shows approximately 10% reduction compared to case-9. The reason is that since the supplies were available at the beginning, they were sent to intermediate nodes close to demand points so the delivery of supplies for the 2nd day can start as soon as the demands appear for that day. Figure 4.17 illustrates the details. 0 10 20 30 40 50 60 70 80 90 100 0 2 4 6 81012141618202224262830323436384042444648 Time (hr) Un sa tis fie d D ema nd (%) 92 Figure 4.16 Percent of unsatisfied demand over time for CASE 9 Figure 4.17 Percent of unsatisfied demand over time for CASE 10 The possible combinations to manage the operations in case-10 are larger than any other case. Consequently, the CPLEX solver went over 2.9 million iterations and it took more than 2 days and 16 hours of CPU time to find the optimal solution. It is clear that a problem with complete set of constraints (if the 0 10 20 30 40 50 60 70 80 90 100 0 2 4 6 81012141618202224262830323436384042444648 Time (hr) Un sa tis fied De ma nd (%) 0 10 20 30 40 50 60 70 80 90 100 0 2 4 6 81012141618202224262830323436384042444648 Time (hr) Un sa tis fie d D em and (%) 93 equity constraints were to be added to the problem) with 2-days of operations, cannot be solved by the commercial solver in a meaningful time period. 4.4. Summary of the preliminary numerical experiments The numerical analysis in this chapter was designed to test the proposed formulation and evaluate the properties of the optimization problem. Ten different case studies were generated based on the same structure of an imaginary hurricane scenario to analyze the effects of the different parameters. In general, the proposed modeling framework produced reasonable outcome. It was able to provide the level of details required in the disaster response logistics at the operational level. For simple cases and small-size problems, the commercial solver was able to find the optimal solutions, however, when the difficult constraints such as equity constraints were added or when the time horizon was extended from 1-day to 2-days, CPLEX commercial solver was unable to deliver good results. It is concluded that better solution algorithms or heuristics are needed to address the larger problem instances or real world size problems. 94 Chapter 5: Solution Approaches In this chapter, first some solution approaches for general integer programming from previous studies in the literature are reviewed in section 5.1. Then in section 5.2, the solution approaches that were specifically used in emergency logistics literature are reviewed in more details. After literature review, in section 5.3 a number of solution techniques are proposed for the mathematical model presented in chapter 3. Two sets of algorithms are proposed to solve the different parts of the problem. In section 5.4 solution algorithms are proposed to solve the hierarchical location finding problem. And in section 5.5, some heuristic algorithms are proposed to solve the general integer vehicle routing problem. In section 5.6, more numerical analysis is performed to further evaluate the robustness of the proposed algorithms. Finally, section 5.7 summarizes the developments in this chapter. 5.1 General Solution Approaches for Integer Programs In General, integer programming problems are very difficult to solve. Over the years, different researchers have proposed several very different solution algorithms. Today, the question is how to select the best approach among the list of available general approaches. Algorithm selection has become an art as some algorithms work better on some specific problem instances. A brief discussion of algorithms is presented in this subsection, attempting to expose readers to their characteristics. More detailed review of integer and combinatorial optimization 95 algorithms can be found in the integer programming literature (e.g. Nemhauser and Wolsey (1999)) Historically, linear programming (LP) has been the base for integer programming (IP) solution approaches. LP was invented in the late 1940's. Those examining LP relatively quickly came to the realization that it would be desirable to solve problems which had some integer variables (Dantzig, 1960). This led to development of algorithms for the solution of pure IP problems. The first algorithms were cutting plane algorithms as developed by Dantzig, Fulkerson and Johnson (1954) and Gomory (1963). Land and Doig (1960) subsequently introduced the branch and bound algorithm. More recently, implicit enumeration (Balas 1965), decomposition (Benders 1962), Lagrangian relaxation (Geoffrion, 1974) and heuristic approaches have been used to solve various integer programs. McCarl and Spreen (1997) suggested the following classification of general algorithms for integer programming problems: 5.1.1 Cutting Planes The first formal IP algorithms involved the concept of cutting planes. Cutting planes iteratively remove parts of the feasible region without removing integer solution points. The basic idea behind a cutting plane is that the optimal integer point is close to the optimal LP solution, but does not fall at the constraint intersection so additional constraints need to be imposed. Consequently, constraints are added to force the non-integer LP solution to be infeasible without eliminating any integer solutions. This is done by adding a constraint forcing the non-basic variables to be greater than a small nonzero value. The simplest form of 96 a cutting plane would be to require the sum of the non-basic variables to be greater than or equal to the fractional part of one of the variables. The cutting plane algorithms continually add such constraints until an integer solution is obtained. Methods for developing cuts appear in Gomory (1963) in more details. Several points need to be made about cutting plane approaches. First, many cuts may be required to obtain an integer solution. For example, Beale (1977) reports that a large number of cuts is often required (in fact often more cuts are required than can be computationally afforded). Second, the first integer solution found is the optimal solution. This solution is discovered after only enough cuts have been added to yield an integer solution. Consequently, if the solution algorithm runs out of time or space the modeler is left without an acceptable solution (this is often the case). Third, given comparative performance with other algorithms, cutting plane approaches have faded in popularity (Beale,1977). 5.1.2 Branch and Bound The second solution approach developed was the branch and bound algorithm. Branch and bound, originally introduced by Land and Doig (1960), pursues a divide-and-conquer strategy. The algorithm starts with a LP solution and also imposes constraints to force the LP solution to become an integer solution similar to cutting planes. However, branch and bound constraints are upper and lower bounds on variables. The branch and bound solution procedure generates two problems (branches) after each LP solution. Each problem excludes the unwanted non- 97 integer solution, forming an increasingly more tightly constrained LP problem. There are several decisions required. One must both decide which variable to branch on and which problem to solve (branch to follow). When one solves a particular problem, one may find an integer solution. However, one cannot be sure it is optimal until all problems have been examined. Problems can be examined implicitly or explicitly. Maximization problems will exhibit declining objective function values whenever additional constraints are added. Consequently, given a feasible integer solution has been found, then any solution, integer or not, with a smaller objective function value cannot be optimal, nor can further branching on any problem below it yield a better solution than the incumbent (since the objective function will only decline). Thus, the best integer solution found at any stage of the algorithm provides a bound limiting the problems (branches) to be searched. The bound is continually updated as better integer solutions are found. The problems generated at each stage differ from their parent problem only by the bounds on the integer variables. Thus, a LP algorithm that can handle bound changes can easily carry out the branch and bound calculations. The branch and bound approach is the most commonly used general purpose IP solution algorithm and it is implemented in many commercial solvers. However, its use can be expensive. The algorithm does yield intermediate solutions which are usable although not optimal. Often the branch and bound algorithm will come up with near optimal solutions quickly but will then spend a 98 lot of time verifying optimality. Shadow prices from the algorithm can be misleading since they include shadow prices for the bounding constraints. A specialized form of the branch and bound algorithm for zero-one programming was developed by Balas (1965). This algorithm is called implicit enumeration. 5.1.3 Lagrangian Relaxation Lagrangian relaxation (Geoffrion (1974), Fisher (1981)) is another area of IP algorithmic development. Lagrangian relaxation refers to a procedure in which some of the constraints are relaxed into the objective function using an approach motivated by Lagrangian multipliers. The basic Lagrangian relaxation problem for the mixed integer program involves discovering a set of Lagrange multipliers for some constraints and relaxing that set of constraints into the objective function. The main idea is to remove difficult constraints from the problem so the integer programs are much easier to solve. IP problems with structures like that of the transportation problem can be directly solved with LP. The trick then is to choose the right constraints to relax and to develop values for the Lagrangian multipliers leading to the appropriate solution. Lagrangian relaxation has been mainly used in two settings: 1) to improve the performance of bounds on solutions; and 2) to develop solutions which can be adjusted directly or through heuristics so they are feasible in the overall problem (Fisher (1981)). An important Lagrangian relaxation result is that the relaxed problem provides an upper bound on the solution to the unrelaxed problem at any stage. Lagrangian relaxation has been heavily used in branch and bound 99 algorithms to derive upper bounds for a problem to see whether further branching down on that branch is worthwhile. 5.1.4 Benders Decomposition Benders decomposition is another algorithm to solve integer programs. This algorithm solves mixed integer programs via structural exploitation. Benders (1962) developed the procedure which decomposes a mixed integer problem into two problems; an integer master problem and a linear subproblem. Then these problems are solved iteratively. Consider the following decomposable mixed IP problem: Maximize FX + CZ s.t. GX ? b1 HX + AZ ? b2 DZ ? b3 X is integer, Z ? 0 Assuming X* is a feasible set of points for integer variables X, then the subproblem for any given X* would be: Maximize CZ s.t. AZ ? b2 - HX* (?) DZ ? b3 (?) Z ? 0 Solution to this subproblem yields the dual variables in parentheses. In turn a "master" problem is formed as follows: Maximize FX + Q X, ?, ?, Q s.t. Q ? ?i (b2 ? HX) + ?i b3 for i = 1,2,3?, p GX ? b1 X is integer , Q is unrestricted 100 This problem contains the dual information from above and generates a new X value. The constraint involving Q gives a prediction of the subproblem objective function arising from the dual variables from the ith previous guess at X. In turn, this problem produces a new and better guess at X. Each iteration adds a constraint to the master problem. The objective function consists of FX + Q, where Q is an approximation of CZ. The master problem objective function therefore constitutes a monotonically non-increasing upper bound as the iterations proceed. The subproblem objective function (CZ) at any iteration plus FX can be regarded as a lower bound. The lower bound does not increase monotonically. However, by choosing the larger of the current candidate lower bound and the incumbent lower bound, a monotonic non-decreasing sequence of bounds is formed. The upper and lower bounds then give a monotonically decreasing gap between the bounds. Benders decomposition convergence occurs when the difference between the bounds is driven to zero. When the problem is stopped with a tolerance, the objective function will be within the tolerance, but there is no relationship giving distance between the variable solutions found and the true optimal solutions for the variables. Convergence will occur in a practical setting only if for every X a relevant set of dual variables is returned. This will only be the case if the subproblem is bounded and has a feasible solution for each X that the master problem yields. This may not be generally true. Also the boundedness and feasibility of the subproblem says nothing about the rate of convergence. The real art of utilizing Benders decomposition involves the recognition of appropriate problems and/or 101 problem structures which will converge rapidly. The procedure can work very poorly for certain structures (Sherali 1981). In general: 1. The decomposition method does not work well when the X variables chosen by the master problem do not yield a feasible subproblem. Thus, the more accurately the constraints in the master problem portray the conditions of the subproblem, the faster will be convergence. 2. The tighter (more constrained) the feasible region of the master problem the better. 3. When possible, constraints should be entered in the master problem precluding feasible yet unrealistic (suboptimal) solutions to the overall problem. The most common reason to use Benders method is to decompose large mixed integer problem into a small, difficult master problem and a larger simple linear program. This allows the solution of the problem by two pieces of software which individually would not be adequate for the overall problem. It should be noted that in Benders decomposition method, the master problem is still an integer program that might be very difficult to solve. 5.1.5 Heuristics Many IP problems are combinatorial and difficult to solve by nature. In fact, the study of NP complete problems (Papadimitrou and Steiglitz (1982)) has shown extreme computational complexity for problems such as the traveling salesman problem. Such computational difficulties have led to a large number of 102 heuristics. These heuristics are used when: a) the quality of the data does not merit the generation of exact optimal solutions; b) a simplified model has been used, and/or c) when a reliable exact method is not available, computationally attractive, and/or affordable. Arguments for heuristics are also presented regarding improving the performance of an optimizer where a heuristic may be used to save time in a branch and bound code, or if the problem is repeatedly solved. Many IP heuristics have been developed, some of which are specific to particular types of problems. For example, there have been a number of traveling salesman problem heuristics as reviewed in Golden et al (1980). Zanakis and Evans (1981) provide a general review of heuristics. Generally, heuristics perform well on special types of problems, quite often coming up with errors of smaller than two percent (McCarl and Spreen (1997)). Zanakis and Evans (1981) provide discussions of selections of heuristics vis-a-vis one another and optimizing methods. 5.1.6 Structural Exploitation Past experiences on IP have indicated that general-purpose IP algorithms do not work satisfactorily for all IP problems. Recently, the most promising developments have involved structural exploitation, where the particular structure of a problem has been used in the development of the solution algorithm. Benders decomposition and Lagrangian relaxation are two examples of structural exploitation. Some problem reformulation approaches and also specialized branch and bound algorithms adapted to particular problems are examples of structural 103 exploitation. The main mechanisms for structural exploitation are to develop an algorithm especially tuned to a particular problem or, more generally, to transform a problem into a simpler problem to solve. The application of such algorithms has sometimes led to spectacular results, with problems with thousands of variables being solved in seconds of computer time (McCarl and Spreen (1997)). Unfortunately, none of the available algorithms have been shown to perform satisfactorily for all IP problems. However, certain types of algorithms are good at solving certain types of problems and a number of efforts have concentrated on algorithmic development for specially structured IP problems. The following section reviews some of approaches used in emergency logistics literature. 5.2 Solution approaches used in emergency logistics literature Chapter 2 provided an extensive review of previous research in the emergency logistics literature. From the number of researches discussed in chapter 2 only four publications are found to have a mathematical model that are partially similar to the mathematical model proposed in this research. In the following paragraphs the solution approaches used in these four lead publications are reviewed. Haghani and Oh (1996) proposed a formulation and solution of a multi- commodity, multi-modal network flow model for disaster relief operations. Their model can determine detailed routing and scheduling plans for multiple transportation modes carrying various relief commodities from multiple supply 104 points to demand points in the disaster area. They formulated the multi-depot mixed pickup and delivery vehicle routing problem with time windows as a special network flow problem over a time-space network. The objective was minimizing the sum of the vehicular flow costs, commodity flow costs, supply/demand storage costs and inter-modal transfer costs over all time periods. Structurally, their model was composed of two network flow problems; one with only real-valued variables and the other with integer variables were connected with a set of capacity constraints called linkage constraints. They developed two heuristic solution algorithms; the first one was a Lagrangian relaxation approach, and the second was an iterative fix-and-run process. The first solution algorithm decomposes the model into two subproblems based on the relaxation of linkage constraints. Lagrangian relaxation is used with penalty for shortage of capacity for linkage constraints. The algorithm was iteratively applied until two subproblems converge. The second solution algorithm was an ad hoc method that fixed integer variables gradually. First all integer variables were relaxed and LP relaxation is solved. Then based on the LP solution, the values of some of the integer variables were fixed to an integer value and the LP was solved again. This process was repeated iteratively until all integer variables are fixed to integer values. Haghani and Oh (1996) solved several instances of numerical problems with both algorithms. For smaller size problems, they showed both algorithms were successful in solving integer problem instances much faster than commercial solvers. They also showed for larger problem instances that the commercial solver 105 was unable to find the optimal solution; both algorithms were able to find close to optimal solution in relatively short CPU times. Comparing the two algorithms, they concluded that the proposed fix-and-run algorithm outperforms the Lagrangian relaxation algorithm both in CPU time and final solution quality. Barbarosoglu and Arda (2004) developed a two-stage stochastic programming model for transportation planning in disaster response. Their study expanded on the deterministic multi-commodity, multi-modal network flow problem of Haghani and Oh (1996) by including uncertainties in supply, route capacities, and demand requirements. The authors designed 8 earthquake scenarios to test their approach on real-world problem instances. Their model is a planning model that does not deal with details required at strategic or operational levels. The model does not address facility location problem or vehicle routing problem. To solve numerical examples, Barbarosoglu and Arda (2004) in the first stage generate random scenarios for supply, demand, and available capacity. In the second stage they used the commercial solver GAMS to solve the resulted network flow problem to minimize the cost. They did not propose any special solution algorithms but used GAMS software to solve the numerical studies. Ozdamar et al. (2004) addressed an emergency logistics problem for distributing multiple commodities from a number of supply centers to distribution centers near the affected area. They formulated a multi-period multi-commodity network flow model to determine pickup and delivery schedules for vehicles as well as the quantities of loads delivered on these routes, with the objective of 106 minimizing the amount of unsatisfied demand over time. The structure of the proposed formulation enabled them to regenerate plans based on changing demand, supply quantities, and fleet size. They developed an iterative Lagrangian relaxation algorithm and a greedy heuristic to solve the problem. The Lagrangian relaxation approach used in Ozdamar et al. (2004) was similar to the one previously discussed in Haghani and Oh (1996) with the only change that Ozdamar et al. (2004) used commercial solver GAMS to solve the linear relaxations. The proposed greedy algorithm solves the network flow problem without considering vehicles to find the best routes for the flow of commodities. Then the algorithm assigns the vehicles to the first available shipment so to minimize the shipment delay. If the vehicles are not available immediately, the shipment is postponed till the earliest available vehicle arrives. The greedy approach is myopic in the sense that the vehicles are independently assigned to the first available job instead of considering the other combinations that might be more rewarding. Comparing the Lagrangian relaxation algorithm and the greedy algorithm in Ozdamar et al. (2004), it was concluded that the greedy algorithm performs faster than Lagrangian relaxation algorithm. However, the greedy algorithm usually resulted larger gaps with global optimal compared to the Lagrangian relaxation. Greedy algorithm did not perform well especially when the capacity was tight that is the usual case in disaster response operations. Finally, Yi and Ozdamar (2007) proposed a model that integrated the supply delivery with evacuation of wounded people in disaster response activities. 107 They considered establishment of temporary emergency facilities in disaster area to serve the medical needs of victims immediately after disaster. They used the capacity of vehicles to move wounded people as well as relief commodities. Their model considered vehicle routing problem in conjunction with facility location problem. The proposed model is a mixed integer multi-commodity network flow model that treats vehicles as integer commodity flows rather than binary variables. Their numerical experiment considered a potential earthquake scenario for the city of Istanbul in Turkey. The numerical problem had 20 nodes, 3 transportation modes, 2 relief commodities and modeled for 8 time periods. They used commercial solver CPLEX 7.5 to solve the IP model. They did not propose any new solution algorithm to the problem however they offered an algorithm to find the itinerary of vehicles from the optimal solution output of CPLEX integer programming solver. They reported that post processing algorithm was pseudo- polynomial in terms of the number of vehicles utilized. Yi and Ozdamar (2007) took the network flow vehicle routing (where vehicles are treated as general integer-valued commodities) and compared it with classic 0-1 vehicle routing. They showed that the general integer formulation is more compact and it is much more efficient for solving. They experienced CPU times ?in seconds? for general integer VRP versus ?in minutes? for classic binary VRP. However in general integer VRP, post processing was needed to extract detailed vehicle routing and pickup or delivery schedules. 108 To summarize, it is shown that in previous publications only a few mathematical models can be found which have relatively similar structures to the model proposed in this research. In those publications, three solution approaches are proposed and tested; Lagrangian elaxation, fix-and-run euristic, and greedy heuristic algorithm. Lagrangian relaxation is successful in proving a bound but it was shown to be the most time consuming algorithm. Greedy heuristic algorithm was shown to be faster compared to Lagrangian relaxation algorithm. However, it lacked in the quality of final optimal solution and resulted in large optimality gaps especially when transportation capacity was limited. Fix-and-run heuristic outperformed Lagrangian relaxation in both categories of speed and solution quality. Fix-and-run heuristic compared to Lagrangian relaxation found the final solution in less CPU time and resulted in smaller optimality gap. 5.3 Solution Techniques for Proposed Mathematical Model The mathematical model proposed in chapter 3 is a complex integrated model. Such an integrated model provides the opportunity for a centralized operation plan that can eliminate delays and assign the limited resources to the best possible use. However, the model is a large-scale mixed general integer programming model and solving such a comprehensive mathematical model is a big challenge. As it is shown in preliminary numerical experiments in chapter 4, the commercial solver was unable to find the optimal solution in a reasonable time. Based on the analysis of the solution techniques for similar models in the literature, it is concluded that exact solution algorithms will not be able to efficiently solve the proposed model. Consequently, the best approach might be 109 designing fast heuristic algorithms that can find near optimal solutions in relatively short computation times. On the other hand, since this model is more complicated than all the previous works in the literature, it would be favorable to structurally decompose this problem to some smaller or easier problems. This model integrates commodity flow problem which is a linear multi- commodity network flow problem with multi-echelon facility location problem which is a binary mixed integer program, and multimodal vehicle routing problem which is a large-scale general integer-valued network flow problem. The ?Idea? is to decompose the problem into smaller or easier problems while taking advantage of the special structures that already exist. The multi-commodity network flow problem is a linear program. LP models are considered easy-to-solve since efficient solution algorithms and commercial solvers exist that can quickly solve large-scale linear programs. The difficult parts are the two integer programming subproblems. In the following sections, a number of heuristic algorithms are proposed to solve the integer programming part of mathematical model. First in section 5.4, four heuristics are proposed to solve the hierarchical location finding problem. Then in section 5.5, four new heuristic algorithms are proposed to solve the general integer vehicle routing problem. 5.4 Algorithms for solving the location problem As discussed earlier, the mathematical formulation presented in chapter 3 is composed of three subproblems. The linear commodity flow subproblem is considered easy and can be solved in conjunction to the facility location problem. 110 On the other hand, the general integer vehicle routing subproblem is a large-scale mixed integer program itself which is considered very difficult to solve. This problem is not mathematically decomposable and it is important to keep the interrelations between the three subproblems. To do so, it is suggested to first relax the integrality condition of vehicle routing subproblem and try to solve the location problem. When the optimal locations are known, it would be much easier to solve the vehicle routing problem. Considering relaxed VRP problem inside the location finding problem is a big advantage because it is easier to solve meanwhile it still reflects the effects of the VRP and available transportation capacity on the location finding problem. The mathematical formulation of this location problem can be obtained by only relaxing the mijtY variables (general integer variables related to vehicle routing problem) in the original model presented in chapter 3. In the following subsections, four solution approaches are considered to solve the location finding problem in the form explained in the last paragraph. 5.4.1 Explicit Enumeration The candidate sites for temporary facility locations are chosen prior to emergency response. Consequently, the number of potential sites is known and the number of possible combinations for facility locations is a finite number. The simplest conceivable optimization approach is explicit enumeration. It is possible to generate all possible solutions, evaluate each of them, and keep the best. To test the applicability of explicit enumeration, let?s use the numerical example introduced in Chapter 4: 111 Combinations for selecting 2 MOB out of 4 candidates: 6!2!2 !4 = Combinations for selecting 2 FOSA out of 4 candidates: 6!2!2 !4 = Combinations for selecting 4 SSA out of 10 candidates: 210!6!4 !10 = Total number of combinations is equal to 756021066 =?? . For any given locations, the remaining problem is a linear program that has a network structure. Linear network problems are considered easy to solve since good algorithms and efficient commercial solvers are developed to solve that problem. For instance, for linear relaxation of the numerical experiment introduced in chapter 4 with given locations, CPLEX solver was able to solve the problem in around 7 seconds on average. If it is required to enumerate all combinations, the total CPU time is equal to hours7.14sec52920sec77560 ==? . It can be concluded that since it is easy to solve the problem after locations are given, it is still possible to explicitly enumerate all combinations and find the final optimal solution. It might not be wise to solve for every single combination, however, it indicates the level of difficulty of the IP problem and provides a benchmark for development and comparison of other solution algorithms. Some other heuristic algorithms are introduced in the following subsections. 5.4.2 Branch and bound - Hierarchical decomposition Branch and bound algorithm is widely used to solve integer programs. It is especially successful when the integer variables are 0-1 binary variables as it is the case in location finding problems. Good algorithms and efficient commercial solvers are developed that use the branch and bound technique. ILOG CPLEX 112 solver is a commercial solver that can apply branch and bound to solve binary mixed integer programs. The proposed mathematical model contains three levels of temporary facilities. Mobilization centers (MOB) are at the top. Federal operational staging areas (FOSA) are the intermediate level facilities that receive commodities from MOB. Then there is state staging areas (SSA) that receive commodities from FOSAs. It is possible to use branch and bound to solve all three levels simultaneously. However, it is possible to hierarchically decompose the facility location problems and solve them consecutively. Three decomposition approaches are proposed and tested: 1. Top-to-Bottom: Decompose the problem into federal level facilities and state level facilities. Assume all state level facilities are open (i.e. SSAiLoc i ??= 1 ). Solve the integer program to find the optimal locations for federal level facilities. Fix the solution for top level to its optimal values and solve the integer program for the state level facilities. 2. Bottom-to-Top: Decompose the problem into federal level facilities and state level facilities. Assume all federal level facilities are open (i.e. MOBFOSAiLoci ???= 1 ). Solve the integer program to find the optimal locations for state level facilities. Fix the solution for bottom level and solve the integer program to find the optimal state level facilities. 3. Tier-by-Tier: First solve the integer program to find the optimal locations for MOB level facilities assuming all other facilities are open. Then fix the 113 optimal MOB, assume all SSA are open and solve IP for FOSA facilities. Finally, fix optimal MOB and FOSA then solve for SSA. Table 5.1 shows the results of applying the abovementioned approaches to the numerical problem in chapter 4. Comparing the total CPU times, it can be seen that Tier-by-Tier decomposition resulted in the least computation time. It was able to reduce the CPU time from 379 seconds when all tiers are considered together, to about 203 seconds (a reduction of about 46%). The Top-to-Bottom approach also gives good results with a total of 215 seconds computation time (43% reduction). On the other hand, it seems that for the current example, Bottom-to-Top approach did not provide favorable results. Mainly, when all federal level facilities are forced to be open, it generates an unnecessarily large number of combinations. Exploring all those combinations result in higher than usual computation times in Bottom-to-Top approach. Table 5.1 Branch and Bound and Hierarchical Decomposition Case Solution Time Final Obj (E+7) Iterations Total Time (S) Solve for All Location Tiers 378.69 3.83595 204402 378.69 All State Level = 1, Solve for FED level 191.91 3.83595 181589 Given FED level, SOLVE for State 22.92 3.83595 43781 ALL FED level = 1, Solve for State 819.23 3.77795 559223 Given State, Solve for FED 138.97 3.83595 106213 Solve for MOB, Rest = 1 151.03 3.82113 139943 Given MOB, Solve for FOSA, SSA =1 28.66 3.83595 59960 Given MOB & FOSA, Solve for SSA 22.94 3.83595 43781 214.83 958.2 202.63 Solution times for solver CPLEX 11.0 on Dell desktop with 3GHz CPU and 4GB RAM It is important to mention that all three proposed approaches provided the same optimal locations. Although it is not a proof, it is a very favorable property 114 to have a number of heuristic algorithms that find the exact solution. The design of proposed hierarchical decompositions allows the heuristics to find the exact optimal solution by not cutting the feasible region. For example in Tier-by-Tier approach, when solving for top tier (MOB level), all other lower level facilities are forced to be open regardless of the limitation on the maximum number of open facilities in lower levels. This provides the chance to find the optimal locations for the tier in hand because all lower levels facilities are at their best theoretical combination. 5.4.3 Highest Capacity Ratio Solving linear relaxation of integer programming problems and analyzing the results can reveal very valuable insights. The idea in this heuristic is to use the linear relaxation to find the facility or facilities that are most important for the performance of the system. Returning to capacity constraints in the mathematical formulation in chapter 3, the following equation enforces the sum of all flows leaving facility i, to be less than the loading capacity of facility i if it is selected to be open; or to be zero otherwise: tmiLocLcapX imit c j cm ijt ,,????? Repeated (3.12) If the binary integer variable Loci is relaxed to take any real number between 0 and 1, it can show the capacity ratio that is used in facility i. The facilities with higher capacity ratios are more favorable because they handle the most flow and their existence is more important to the entire response operations. 115 The steps of the Highest Capacity Ratio (HCR) Algorithm are: Step 1- Relax the integrality condition for all temporary facility variables Step 2- Add 10 ?? iLoc for all relaxed binary variables Step 3- Solve the linear relaxation problem and obtain optimal values for all Loci variables Step 4- For each facility type; sort the Loci variables in descending order Step 5- For each facility type; select the facilities with highest capacity ratio from top of the list until the maximum number allowed is reached By following these five steps, one can find the selected facilities in a single snapshot. However, it can be argued that selecting a facility may affect the selection of others. So it might be beneficial to select the one facility with the highest ratio, solve the linear relaxation again, and repeat until maximum number of each facility is selected. The steps of Iterative Highest Capacity Ratio (IHCR) algorithm are: Step 1- Relax the integrality condition for all temporary facility variables Step 2- Add 10 ?? iLoc for all relaxed binary variables Step 3- Solve the linear relaxation problem and obtain optimal values for all Loci variables Step 4- Find the facility i with the highest Loci value Step 5- If the maximum number of facilities are not reached, select facility i, add 1=iLoc to the formulation and go to Step 3. Otherwise stop. 116 To test HCR algorithm, the LP relaxation of the numerical experiment formerly introduced in chapter 3 is solved again. Table 5.2 shows the values for relaxed Loci variables. The constraints for the maximum number of facilities required the selection of 2 MOB out of 4, 2 FOSA out of 4 and 4 SSA out of 10 potential SSA nodes. Solving the linear relaxation with additional 10 ?? iLoc constraints only took about 32 seconds. The resulted node selection shown in Table 5.2 is obtained with using the single snapshot HCR algorithm. It is worth mentioning that for this example, the HCR algorithm was able to find the exact optimal solution and did so incredibly faster than branch and bound method (32 seconds versus 379 seconds). Table 5.2 Values of Loci variables for HCR heuristic algorithm Facility Type MOB FOSA SSA (node number) Loci Value (4) 0.6807 (8) 1.0 (12) 1.0 (17) 0.6174 (5) 1.0 (9) 0 (13) 0.5357 (18) 0.3131 (6) 0 (10) 0.8532 (14) 0.4164 (19) 0.8359 (7) 0.3193 (11) 0.1468 (15) 0 (20) 0.2054 (16) 0 (21) 0.0762 Selected Nodes 4 , 5 8 , 10 12 , 13 , 17 , 19 5.4.4 Static Network-Location Considering a time varying structure and a time-space network is essential to capture the details of emergency response logistics at the operational level. However, it expands the size of the formulation drastically and makes the problem extremely difficult to solve. The idea for this heuristic is to build a static version of the formulation that can be solved much easier and faster. It should still 117 consider the special structure of the network and account for supplies, demands, and facility capacities; but manage to aggregate over the time dimension in order to generate a smaller formulation. To do so, the following mathematical formulation is proposed: ??? m c ji cm ij m ij XtMin , . (5.1) cUiSupXX ci m j cm ji m j cm ij ,???? ???? (5.2) cWiXX m j cm ji m j cm ij ,0 ??=? ???? (5.3) cViDemXX ci m j cm ji m j cm ij ,??=? ???? (5.4) mUiLcapX mi c j cm ij ,????? (5.5) mVUiUcapX mi c j cm ji ,+????? (5.6) mWiLocLcapX imi c j cm ij ,. ????? (5.7) mWiLocUcapX imi c j cm ji ,. ????? (5.8) maxLocLoc i i ?? (5.9) ( ) WiLoci ??? 1,0 and mcjiX cmij ,,,0 ?? The notations are similar to the original problem that is previously defined in section 3.5 with the exception that time index t is dropped from all variables and parameters. As a result, all variables and parameters are static and defined as the aggregate value of the original variables over all time periods. For example, ciSup 118 and ciDem are the aggregate supply and aggregate demand of commodity c in node i, over the entire planning horizon. Decision variable cmijX is the aggregate amount of commodity c that is shipped from node i to node j with transportation mode m, over the entire planning horizon. In this new formulation the details of unsatisfied demand over time is not available. Consequently, the objective function (5.1) is chosen to minimize the total travel time by all commodities. Equation (5.2) and (5.4) enforce the supply and demand constraints for each node and each commodity. Equation (5.3) imposes the conservation of flow at intermediate nodes. Loading and unloading capacity constraints are defined in equations (5.5) and (5.6) for the permanent facilities. Similar constraints for temporary facilities are required by equations (5.7) and (5.8). Finally, equation (5.9) enforces the maximum number of open facilities for each facility type. In the proposed static formulation, vehicle routing constraints are dropped from the formulation. It is equivalent to assume that ample transportation capacity is available or the initial distribution of vehicles is done in such a way that does not affect the choice of temporary facilities. Also, time-space structure is removed from the original model. It can be explained if the variations of supplies, demands, and capacities over the planning horizon are not very drastic. No link capacity is imposed in this formulation; however capacity limitations are reflected in loading and unloading capacities for each facility. It should be noticed that the static formulation is still an integer programming model. However, it is of much lower size and complexity compared 119 to the original formulation while still reflecting the structure and important properties of the original model. Similar to previous heuristics, static network location problem (SNLP) heuristic is also tested with the numerical example of chapter 4. CPLEX solver version 11.0 is used to solve the problem on a Dell desktop computer with 3 GHz CPU and 4GB of RAM. After presolve modifications, reduced MIP had 130 rows, 491 columns, and 1713 nonzeros. It took less than 1 second to solve the modified problem which is extremely faster than the previous heuristics. However, optimal locations obtained from this formulation do not exactly match with the optimal locations of the original IP problem. Using the locations suggested by SNLP results in 2.5% higher objective function value compared to the case that exact optimal locations are used. To summarize, four heuristic approaches are proposed to solve the location finding problem. Computation times vary greatly across the algorithms ranging from 14 hours to less than 1 second. Firstly, explicit enumeration showed that even though LP solution when locations are given takes only 7 seconds, the large number of possible combinations makes it very difficult to explore all the combinations. Secondly, hierarchical decomposition approach suggested that it is beneficial to choose it over the general branch and bound (46% faster). Among the three suggested hierarchical decompositions, Tier-by-Tier decomposition was the fastest. Thirdly, highest capacity ratio heuristic was the fastest among all other heuristics that could still find the exact optimal solution. And finally, SNLP proposed a new formulation that is very efficient and can be solved to find the 120 locations for the original problem. SNLP was the fastest algorithm but the resulted locations were different from those of the exact optimal solution. 5.5 Algorithms for solving Vehicle Routing Problem In section 5.2 the relevant literature that suggested solution methods was summarized. Mainly, three heuristic approaches were proposed to solve the general integer vehicle routing problem: Lagrangian relaxation, Fix-and-Run algorithm, and a greedy algorithm. Using the numerical results, it was also concluded that the Fix-and-Run algorithm proposed by Haghani and Oh (1996) had the best performance. It was the fastest algorithm and it had the least optimality gap compared to the other algorithms. In the following subsections, four heuristic algorithms are proposed to solve the general integer vehicle routing problem. The general idea is adopted from the successful experience of Fix-and-Run heuristic algorithm suggested by Haghani and Oh (1996). The main steps of the proposed algorithms are: 1. The mixed integer linear problem is solved with the relaxation of integer variables. 2. The values of some integer variables are fixed in an orderly manner and the problem is solved again with the relaxation of the remaining integer variables iteratively. 3. When all integer variables are fixed, the process is terminated. 5.5.1 T-Counter Heuristic The steps of T-Counter algorithm are as following: 121 Step 1- Relax all general integer variables and solve the relaxed LP. Set t=0 Step 2- Check all mijtY variables for current time period t. If all mijtY variables are integer, then if t = tlast, terminate. Otherwise, set t = t + 1 and restart Step 2. Step 3- For current time period t, fix all mijtY variables to the closest integer number Step 4- Create a new problem by adding ( mijtY = the fixed value from step 3) constraints to the problem Step 5- Relax the rest of the integer variables and solve the new LP problem Step 6- Set t = t + 1 and go to Step 2 In this algorithm, starting from the first time period, Y variables are fixed iteratively and in a chronological order. If the flow of the vehicles through the network is fixed to be integral at time period t, because of the network structure of the problem, it is more likely that the flows at time periods after t, also turn out to be integral. Conservation of flow in a time-space network requires that if the flows that enter a node are integer, then sum of the flows that leave the same node must also be integral. This does not mean that every single flow leaving that node will definitely be integer but it is a necessary condition. If the planning horizon of the problem is consisted of T time periods, then in worst case the algorithm will go through only T iterations. It is the worst case scenario and not the average case because during an iteration if all Y variables are already integer, the algorithm directly proceeds to next t without solving a LP relaxation. This is a very important property to have because this algorithm will 122 stop at most after T iterations. Fast convergence rate is expected from this algorithm. 5.5.2 Origin-Based T-Counter Heuristic In the previous algorithm, in each iteration all the mijtY variables for current time period t are fixed at the same time. That approach reduces the flexibility of the algorithm to reroute the vehicles within one time period which can sometimes cause suboptimal assignments. To remedy this, Origin-based T-Counter heuristic algorithms is proposed. In this algorithm, outgoing flows from only one origin node will be fixed during each iteration. In other words, for current time period t, we start from node i = 1 and fix all outgoing mijtY variables, solve LP relaxation, then fix all flows from node i = 2, and move to the next node until all nodes are fixed. Then the same procedure is followed for the next time periods until the end of the planning horizon. The steps of Origin-based T-Counter algorithm are: 1- Set t = 0 and i=1 2- Relax all general integer variables and solve the relaxed LP. 3- If all relaxed variables are integer, the IP solution is found, Terminate 4- For current t and i, fix all mijtY variables to the closest integer number 5- Create a new problem by adding ( mijtY = the fixed value from step 4) constraints 6- If i < ilast then set i = i + 1 and go to step 2. Otherwise go to the next step 7- If t = tlast terminate otherwise set i = 1 and set t = t +1, go to step 2 123 This algorithm is more general compared to T-Counter algorithm. If the planning horizon of the problem is consisted of T time periods and N is the number of nodes in the network, then at most T ? N iterations are required to solve the problem. Again this is a worst case scenario and in general the algorithm is expected to find the integer solution before going through all T ? N iterations. 5.5.3 Y-List Heuristic In the previous two algorithms, several Y variables are fixed during each iteration. For example in T-Counter algorithm, at first iteration all mijtY variables with t = 0 are fixed simultaneously that can lead to under utilization of the available vehicles. For more clarification assume a hypothetical scenario where there are 4 vehicles available at node i and 3 exactly similar arcs are leaving node i. Solving the linear relaxation of the problem will assign 1.33 vehicles to each path. Applying T-Counter algorithm or even Origin-based T-Counter algorithm to this example rounds down 1.33 and as a results it assignments 1 vehicle to each path and 1 vehicle will remain unused. The idea of Y-List algorithm is to fix this problem by only selecting one Y variable in each iteration. This will allow the LP model to adjust itself and take advantage of any potential vehicles that might be available and are not being used due to rounding down. Returning to our hypothetical scenario, if the 3 arcs are fixed one by one then all available vehicles will be used. The vehicle assignment will be 1, 1, 2 and all 4 vehicles are utilized. 124 To run this algorithm, it is required to have a priority list of all Y variables. When the first LP relaxation is solved, the algorithm needs to select a Y variable among all non-integer Y variables to fix. It is faster to have a pre- populated list of all Y variables and then fix them one by one if they have a non- integer value. The steps of Y-List algorithm are: 1- Populate a sorted list of all mijtY variables 2- Relax all general integer variables and solve the relaxed LP 3- If all mijtY variables are integer, save the solution & terminate the algorithm, otherwise 4- Select the 1st mijtY from the list, Fix it to the closest integer number and remove it from the Y-list 5- Create a new problem by adding ( mijtY = the fixed value from step 4) constraint 6- Go to step 2 Theoretically in the worst case scenario, the algorithm can go through Y iterations. Y is the total number of all mijtY variables and also the size of the Y- List set. In large scale numerical problems, Y can be a very large number. For example in the numerical experiment in chapter 4, thousands of mijtY variables exist. In the worst case scenario the algorithm need to go over thousands of iterations and fix every single Y variable. However, as it will be shown, usually the algorithm does not need to fix every single Y variable before finding an IP 125 solution. In fact due to having a network structure, an IP solution is found very quickly and the algorithm converges relatively fast in typical numerical examples. 5.5.4 Y-List Modal Heuristic In large-scale logistic operations often multiple transportation modes are utilized. From theoretical perspective, each transportation mode can be considered as the flow of a special commodity over the network. Different transportation modes are not competing for share resources and there is no explicit constraint that relates the flow of different modes. Consequently, it can be assumed that each transportation mode is acting somehow independently. It should be mentioned that relief commodities that are carried by each transportation mode can be transferred to another mode inside intermodal terminal but the vehicles of each mode are never interchangeable. For example, if 2 trucks and 2 helicopters enter a node, always the same 2 trucks and 2 helicopters have to leave that node and it can never transform into 3 trucks and 1 helicopter. Taking advantage of this independence among multiple transportation modes is the idea behind Y-List Modal heuristic. Y-List Modal is very similar to previously described Y-List heuristic, however it tries to fix a Y variable from each transportation mode during each iteration. For example, if two transportation modes exist, the algorithm will fix two Y variables in each iteration. Consequently, if M is the number of available transportation modes, the algorithm will fix M variables in each iteration and it can stop after Y / M iterations in the worst case. 126 The steps of Y-List Modal algorithm are: 1- Populate a sorted list of all mijtY variables for each mode m 2- Relax all general integer variables and solve the relaxed LP 3- If all mijtY variables are integer, save the solution & terminate the algorithm, otherwise 4- For each mode m, Select the 1st mijtY from the list, Fix it to the closest integer number and remove it from its Y-list 5- Create a new problem by adding ( mijtY = the fixed value from step 4) constraints 6- Go to step 2 5.5.5 Comparing Performance of the Proposed Algorithms In previous sections, four heuristic algorithms are proposed to solve the general integer vehicle routing problem. In this section, these algorithms are analyzed and their performance is compared. All four algorithms are applied to the similar numerical example that is previously defined in Chapter 4. The facility location problem is solved in advance and the optimal locations of the facilities are assumed to be known at this stage. The mathematical model is generated and initially solved by CPLEX Software. Table 5.3 represents the statistics of the mathematical model and also the optimization results obtained by the commercial solver. It is shown that the problem is a large-scale mixed integer program with a large number of general integer variables. CPLEX version 11.0 is used on a Dell desktop computer with 3 127 GHz Intel CPU and 4 GB of RAM. As it can be seen in the table, the commercial solver got to as close as 0.5 percent gap but it was unable to find the exact solution for the problem even after a long computation time. 0.5 percent optimality gap should be acceptable in many applications; nonetheless it shows the difficulty of solving the MIP problem even with a strong commercial solver on a fast computer. Table 5.3- Model Statistics and Optimization Results from CPLEX Solver Problem Stats Objective nonzeros = 3881 No. of Variables : 110572 [Nonnegative : 48300, Binary : 18, General Integer : 62254 ] No. of Linear Constraints : 36593 [Equality : 11960 , Non-equality : 24633 ] Nonzeros : 372305 [RHS nonzero : 1467] CPLEX Optimization Results Objective Value(E+7) Solution Time(s) GAP (%) Initial LP Bound MIP Best Bound(E+7) Comments 3.8709 81000 0.51 3.8059 3.8511 User-Stopped after 22.5 hrs 3.9121 1041 1.98 3.8059 3.8345 Stopping gap set to 2% Table 5.4 shows the results of solving the same problem using the four heuristic algorithms proposed in this chapter to solve general integer vehicle routing problem. Comparing gap percentiles from the best IP, it is concluded that the proposed algorithms were generally successful. Three of the four proposed heuristic algorithms provided very small optimality gaps of between 1 and 2.5 percent to the best IP solution provided by the commercial software after 22.5 hours. Comparing the solution times is even more impressive. It can be seen that 128 all algorithms found an IP solution and all of them converged in less than about 4 minutes. It is very important to quickly find close to optimal solutions especially in this problem that deals with dynamic emergency response operations. Table 5.4 Numerical results of the proposed VRP Heuristic Algorithms Algorithm Objective Value (E+7) % GAP (Initial LP) %GAP (Best IP) Iterations Solution Time (s) T-Counter 4.2525 10.85 9.85 98 113.3 Origin-Based T-Counter 3.9668 3.41 2.47 3977 247.1 Y-List 3.91615 2.09 1.16 851 89.1 Y-List Modal 3.9300 2.45 1.52 507 73.7 Comparing the four algorithms, the Y-List algorithm is shown to find the best solution quality with the minimum gap. Y-List Modal and Origin-Based T- Counter algorithms also resulted in very good objective functions and small optimality gap. T-Counter algorithm has the largest gap of about 10 percent. It should be reminded that the idea for T-Counter algorithm was adopted from Haghani and Oh (1996) which was the best practice in literature available to this date, to the best of our knowledge. Comparing the solution speed and rate of convergence, it can be seen that all algorithms are quite fast. Y-List Modal was the fastest algorithm with only 73.7 seconds CPU time. Y-List and T-Counter algorithms are in 2nd and 3rd place with relatively close solution times. Y-List Modal produced the longest solution time of about 4 minutes, mainly due to the large number of iterations that was required. It is very important to notice that the number of iterations is not directly 129 related to the solution time, because different iterations take different CPU times. For example, Origin-Based T-Counter goes through about 4000 iterations in about 4 minutes compared to about 100 iterations of T-Counter that takes about 2 minutes. Also, Y-List Modal that recorded the least solution time does not have the least number of iterations. Figure 5.1 shows the convergence rate of the four algorithms. All algorithms initially start from LP relaxed solution which is an infeasible solution for the IP problem. Over time, algorithms try to find integer solutions and reduce this infeasibility. As more and more integer variables are found, the objective function increases. In this way, as soon as an all-integer set of variables are found; the algorithms will stop and report the best solution that is feasible for the IP problem. Figure 5.1 shows a steep slope only for T-Counter algorithm and all other algorithms have a steady and very gradual slope. The main reason is that T- Counter algorithm fixes a large number of integer variables in every iteration that reduces the number of iterations but on the other hand does not permit the LP relaxation to re-adjust and utilize the vehicles that are left behind due to rounding down. All other three algorithms, fix a very small number of variables in every iteration. This allows the LP relaxation to adjust to the fixed values and re-route the commodities and vehicles to take advantage of any remaining transportation capacity. 130 37 38 39 40 41 42 43 44 45 0 60 120 180 240 Ob jec tive Fun ct ion (m illi on s) Elapsed Time (sec) T-Counter Origin-Based T-Counter Y-List Y-List Modal BEST LP BEST IP Figure 5.1 Convergence rate comparisons of the proposed algorithms Table 5.5 summarizes the analysis and comparisons of the four proposed algorithms. Each row shows a criterion and comparatively ranks the four algorithms for those criteria. For example in best solution criteria, Y-List and Y- List Mode are ranked one and two. Comparing the convergence rate, it can be seen that Y-List Modal was the fastest algorithm followed by Y-List algorithm. On the other hand when the least number of iterations are compared, T-Counter is the winner. Also for theoretical worst-case criteria, T-Counter and Origin-Based T-Counter algorithms are ranked 1st and 2nd. As explained earlier, one iteration of each algorithm does not take the same amount of time as one iteration of the other algorithms. Origin-Based T-Counter has the fastest time per iteration followed by Y-List algorithm. 131 Table 5.5 ? Comparative Rankings of the Proposed Heuristic Algorithms Comparative Ranking T-Counter Origin-Based T-Counter Y-list Y-list Modal Best Solution Quality 4th 3rd 1st 2nd Convergence Speed 3rd 4th 2nd 1st Least No. of Iterations 1st 4th 3rd 2nd Best worst- case Scenario 1st 2nd 4th 3rd CPU Time per Iteration 4th 1st 2nd 3rd To summarize, it is shown that all four algorithms are capable of finding good quality solutions in relatively short computational times. Having short computation time is the most important property of the proposed algorithms that makes it possible to apply them in real-world dynamic operations. In Fact, the applicability of proposed mathematical model in chapter 3 could be hardly justified without fast solution algorithms that can adjust and re-optimize in real- time. Comparing the four algorithms, it is concluded that no single algorithm dominates the others in all criteria rankings. When solution quality and convergence speed is more important, Y-List and Y-List Modal are showed to perform better. On the other hand, when good performance under worst-case scenario is important, T-Counter and Origin-Based T-Counter algorithms are shown to have better statistics. It should be noted that all of the four proposed algorithms are heuristic algorithms. Even though they showed very impressive results for the current 132 numerical experiment, there is no proof that they will always have equally good performances for all problem instances. As explained in section 5.1.5, this is in the nature of most heuristic algorithms and is not limited to this study. However, to test the robustness of the proposed algorithms further, more numerical experiments are conducted in the following section. 5.6 Testing robustness of the proposed VRP heuristics In section 5.5, four heuristics solution algorithms were proposed to solve the general integer vehicle routing problem. All four algorithms showed reasonably good result for a particular numerical experiment. In this section, more cases are generated and solved to test the robustness of the proposed algorithms in various conditions. A total of ten random numerical cases are generated and solved by all four proposed heuristics algorithms as well as the CPLEX commercial solver. In all these ten cases, the network structure is similar to FEMA?s structure that was introduced in previous numerical experiments. Also, for all these cases the location problem is solved in advance and the optimal locations of all levels of temporary facilities are known and fixed. Table 5.6 lists the objective function values for these 10 cases. Linear relation of the problem is also solved and the objective function of relaxed LP is reported as a lower bound for comparison. Also in Table 5.6, average of the objective functions for the 10 cases and its standard error is reported. The last row of the table shows the average gap between the final solution of each algorithm and the optimal solution for linear relaxation of the problem. 133 Table 5.6 Objective function values for 10 random cases Objective Function LP relaxation T-Counter Origin-bsd T-Counter Y-list Y-list Modal CASE 1 3.50445 3.79365 3.5454 3.52815 3.54615 CASE 2 3.46545 3.69225 3.5121 3.4965 3.48285 CASE 3 3.5199 3.7290 3.56865 3.5481 3.5418 CASE 4 3.53715 3.84465 3.5979 3.5712 3.56535 CASE 5 3.49095 3.7296 3.5433 3.52395 3.51495 CASE 6 3.48825 3.73365 3.55515 3.51705 3.51255 CASE 7 3.73995 3.98055 3.77055 3.75705 3.7791 CASE 8 3.4392 3.7134 3.4773 3.45825 3.4677 CASE 9 3.50805 3.7758 3.55755 3.5331 3.5388 CASE 10 3.5262 3.76065 3.5766 3.5568 3.5586 Average 3.521955 3.77532 3.57045 3.549015 3.550785 Coef.of.Var 0.0232 0.0223 0.0218 0.0224 0.0242 Avg GAP(%) 0 7.20 1.38 0.77 0.82 Across the 10 random cases, all four algorithms present consistent performances in general. Comparing the average gap percentile, it is shown that Y-list modal had the best overall performance closely followed by Y-list modal. For both of these cases the average gap is less than 1 percent which is a very favorable outcome. Origin-based T-counter algorithm has also acceptable results with only 1.38% optimality gap on average. The largest gap is produced by T- counter algorithms at 7.2 percent on average. 134 Table 5.7 lists the running time of the proposed algorithms for the same 10 random cases. On average, Y-list algorithm has the fastest convergence rate at about 2 minutes. The running time for other 3 algorithms is also in the acceptable range. Highest average running time belongs to Origin-based T-Counter algorithm at about 5 minutes. Table 5.7 Running time comparisons of the algorithms for 10 random cases CPU Time (sec) T-Counter Origin-based T-Counter Y-list Y-list Modal Case 1 191 210 119 159 Case 2 107 197 137 84 Case 3 170 420 155 140 Case 4 284 481 132 112 Case 5 263 270 110 97 Case 6 211 357 111 127 Case 7 192 399 124 307 Case 8 254 266 109 161 Case 9 195 262 92 172 Case 10 243 304 179 162 Average 211 316.6 126.8 152.1 Coef.of.Var 0.2468 0.2971 0.1998 0.4080 5.7 Summary In this chapter, first some solution approaches for general integer programming from previous studies in the literature were reviewed. It was concluded that the current model is very complex and a reliable exact solution method is not available that would be computationally attractive or affordable. 135 Consequently, it is more realistic to develop fast and efficient heuristic algorithms to find near optimal solutions. The solution approaches exclusively used in emergency logistics were also reviewed. Most studies used commercial solvers to solve their model, only three solution approaches were proposed and tested; Lagrangian Relaxation, Fix-and- Run Heuristic, and Greedy Heuristic algorithm. Lagrangian relaxation was successful in proving a bound but it was shown to be the most time consuming algorithm. Greedy Heuristic algorithm was shown to be faster compared to Lagrangian relaxation algorithm but the quality of final optimal solution was not good. Fix-and-Run heuristic outperformed Lagrangian Relaxation in both categories of speed and solution quality. Fix-and-Run heuristic compared to Lagrangian relaxation found the final solution in less CPU time and resulted in smaller optimality gap. To solve the mathematical model in this research, it was structurally decomposed into three sub-problems: multicommodity network flow, location finding with multiple layers, and general integer vehicle routing problem. These 3 problems were solved one-by-one, however the interrelation between these 3 problems were preserved at all times. Multicommodity network flow problem is a linear program and is considered easy since efficient commercial solvers exist that can solve very large LP programs quickly. To solve the multi-layer facility location problem, 4 heuristic methods are proposed. From those, the Branch-and-bound with 136 hierarchical decomposition and the highest capacity ratio were the 2 algorithms that showed better results. To solve the general integer vehicle routing problem four heuristic algorithms were proposed. The algorithms were tested with a large size numerical experiment. All four algorithms were successful in finding a good integer solution. The convergence rates of the proposed algorithms were also much faster than the commercial solver for the same optimality gap. The proposed VRP algorithms were compared to each other. It was concluded that Y-list and Y-list Modal algorithms were better in solution quality and the convergence speed. However, when worst-case scenario is considered, T- counter and Origin-based T-counter algorithms were shown to have better performances. Finally, all four algorithms were used to solve 10 random generated problem instances. It was concluded that the proposed algorithms could find good solutions very quickly. In fact, for most cases less than 2% optimality gap was reached in less than 2 minutes of computation time. 137 Chapter 6: Detailed Analysis of the Mathematical Model In this chapter, in-depth analyses of different aspects of the proposed mathematical model are presented. These analyses are provided in order to better illustrate the capabilities of the model and also examine model?s sensitivity in various circumstances. The analyses in this chapter are divided into two main categories: sensitivity analysis of the structural parameters of the model, and sensitivity analysis of the main input values of the model. In the following subsections, first in section 6.1 a numerical case study is introduced that is being used to perform the analyses in the rest of this chapter. Then in section 6.2, some parameters that affect the structure of the model are introduced and sensitivity analysis is performed to investigate their role. In section 6.3, sensitivity analysis is performed over several input parameters. It is shown that changing some input parameters not only affects the optimization results but it can also largely alter the problem size and solution computation times. Section 6.4 summarizes the overall findings of the sensitivity analysis. 6.1 Introduction of the numerical case study The numerical problem in this chapter is an imaginary scenario where a natural disaster such as a hurricane strikes the southern coast of the United States. It is assumed that two separate regions, one in Mississippi and one in Louisiana, are affected. The network structure of the problem is similar to the case study introduced in chapter 4. One logistics center (LC), one commercial storage site (CSS) and one vendor (VEN) are the three main permanent sources to store and 138 ship the relief items. Three levels of temporary facilities that receive and transfer relief items and vehicles include: Four mobilization centers (MOB), four federal operational staging areas (FOSA) and ten state staging areas (SSA). The demands are concentrated at twenty points of distribution (POD) between two disaster areas. Figure 6.1 shows the disaster area and the locations of the facilities. The other definitions and parameter values are similar to those expressed in chapter 4 unless otherwise is stated in the following subsections. Figure 6.1 Disaster area map and facility locations The computer used for the computational experiments presented in the rest of this chapter is a Dell desktop computer with Intel Xeon 3.0 GHz CPU with 3.5 GB of RAM and Windows XP operation system. To solve the mathematical formulation, ILOG CPLEX 11.0 is used. Microsoft Visual Basic 6 is used to create the formulation and post processing of the optimization results. LC(1) CSS(2) VEN(3) MOB(4)MOB(5)MOB(6) MOB(7) FOSA(8) FOSA(9) FOSA(10)FOSA(11) SSA(12) SSA(13) SSA(14) SSA(15) (22) (23)(24) (25)(26) (27) (28)(29) SSA(16) SSA(17) SSA(18) SSA(19) SSA(21) SSA(20) (30) (31) (32) (33) (34) (35) (36) (37) (38) (39)(40)(41) 139 6.2 Sensitivity Analysis of the Structural Parameters Structural parameters are those parameters in the model that affect the structure of the mathematical formulation. Different values of these parameters can drastically change the size and behavior of the model. For current formulation, some examples of the structural parameters include: 1. Number of commodities C 2. Number of transportation modes M 3. Time-step resolution t In the following subsections, the effects of different values of these parameters are investigated: 6.2.1 Analysis of Number of Commodities In the proposed model, different values of C (number of commodities), fundamentally change the structure of the model and can largely affect the size and the difficulty of solving the model. When C=1, the formulation represents a single-commodity problem but when C > 1 the formulation transforms into a multi-commodity problem that can be more difficult to solve. A multi-commodity problem compared to a single commodity problem requires more data, larger number of decision variables and a larger number of constraints. From decision maker?s perspective, in multi-commodity problems another dimension is added to the problem in order to find the optimal balance between the amounts of several commodities that are being loaded in each shipment. 140 To test the effect of the number of commodities, four cases are considered each with one, two and three and four commodities respectively. In order for all four cases to be relatively comparable, the amounts of total demand and total supply of all commodities at each location are kept the same. However, for any given node the supply and demand is different for each commodity. For example, in two-commodity problem, supply and demand for 1st commodity is assumed to be twice the supply and demand for the 2nd commodity; and so on for three- commodity and four-commodity problems. Using a customized Visual Basic code, the mathematical formulation for these four cases are generated. Table 6.1 reports the problem size for each of these cases. Then each case is solved with CPLEX commercial solver. Since the supply and demand amounts are the same, the objective function values are not that different for all four cases. However, the CPU time to solve each case is different. It takes only 7.11 seconds to solve the single commodity problem. However, for 2, 3, and 4 commodity problems solution times rapidly increase to 71.76, 212.28 and 582.38 seconds respectively. Table 6.1 Problem sizes for different number of commodities Case Description Number of Variables Number of Constraints None-Zero Coefficients Input file size (Kb) Single Commodity 89185 32601 260075 5384 Two Commodities 110570 36585 372297 7355 Three Commodities 134719 40569 484519 9326 Four Commodities 158868 44553 596741 11298 141 Figure 6.2 illustrates the solution time as well as problem sizes for the four cases with different number of commodities. The stopping criteria for optimization is set to be 1% optimality gap. As it can be seen in the graph, when the number of commodities increases, number of variables, constraints and CPU- time also increase. It is interesting to notice that increase in problem size is almost linear; however the CPU-time increases much faster. It can be concluded that larger number of commodities makes the problem exponentially difficult to solve. 0 100 200 300 400 500 600 0 20 40 60 80 100 120 140 160 180 0 1 2 3 4 5 Thou sands Number of Variables Number of Constraints CPU Time se cond Figure 6.2 ? Problem size and solution time versus the number of commodities 6.2.2 Analysis of the Number of Transportation Modes Another example of the structural parameters is M (number of transportation modes) that can have a major affect on the structure of model as well as the difficulty of solving the problem and the behavior of the results. When M=1 only one transportation mode is used to deliver the relief commodities. In a single modal problem, all shipments are transported by only one mode from 142 origins to their destinations. On the other hand, in multimodal systems there is this question of which transportation mode to utilize and how to balance the commodity flows among different modes. Another concern in multimodal environments is considering the intermodal transfer. From application?s perspective it is important to provide suitable facilities and required equipments in order to transfer the relief commodities between transportation modes quickly and efficiently. From modeling perspective it is important to consider the properties of each transportation mode and correctly model the delays during intermodal transfers. Also for intermodal transfer facilities, it is important to consider the loading and unloading capacities based on the availability of the relevant equipments. The mathematical formulation presented in chapter 3 is capable of modeling multiple modes of transportation. It is assumed that main FEMA facilities at federal level have access to multiple modes and also act as intermodal transfer nodes. Equation (3.2) controls the flow of commodities by each mode and also keeps track of commodity transfers between modes. For example mcmitXT ? is equal to the amount of commodity type c in node i which is transferred from mode m to mode m? at time t. Also mmK ? is used as intermodal delay which is equal to the time required to transfer commodities from mode m to mode m? . The numerical example of chapter 4 only considered one mode of transportation which was the ground transportation and only one kind of vehicle which was 53ft trailer truck. To analyze the effect of multimodal operations in this chapter, air transportation mode is added to the problem. For the sake of the 143 numerical example in this chapter, the aircraft of choice is selected to be C130- Hercules cargo planes. C130 has the capacity of about 4500 cft and is assumed to travel at an average speed of 250 mph. Adding a new transportation mode, adds a new layer of network to our time-space framework. Since having air transport facilities at local level or very close to disaster areas might not be possible, it is assumed that only federal level facilities have access to air transportations. In other words, federal level facilities (e.g. LC, CSS, VEN, MOB, and FOSA) are connected through air and ground transportation but state staging areas (SSA) and points of distribution (POD) are only accessible by ground transportation. Consequently, a shipment from logistics center (LC) can be sent with airplane to a federal operational staging area then transferred to ground transportation and then sent to SSA and finally delivered at PODs. This is only an assumption for current numerical example and not a general limitation for the proposed mathematical model. Based on the above description, two numerical cases are formulated and solved. The first case is a single-modal problem with only ground transportation. Supplies, demands, capacities and other parameters are fixed and similar to those in the previous sections. The number of tractor trailers used is 60 which are initially located at the three source nodes (e.g. LC, CSS, VEN) with 20 trucks at each location. In the second case, air-transportation mode is also available between federal-level facilities. Twenty C130 cargo planes with 4500 cft capacity are utilized. Ten planes are initially located at CSS and 10 planes are located at VEN facilities. 144 Table 6.2 represents the problem size for each of these cases. In this example, by adding the second transportation mode, number of variables increases by about 6 percent and number of constraints increases by about 8 percent compared to the single modal case. Both cases are solved by CPLEX commercial solver and Table 6.3 summarizes the optimization results. Comparing the objective functions (sum of all unsatisfied demand over demand nodes, commodities and time periods) it is evident that by introducing a new transportation mode, the objective function is improved. It is expected because new transportation mode increases the transportation capacity which results in faster delivery of relief commodities. Table 6.2 Problem sizes for different number of transportation modes Case Description Number of Variables Number of Constraints None-Zero Coefficients Input file size (Kb) Single Modal 104641 33902 349785 6966 Multimodal 110570 36585 372297 7355 The improvement in objective function value is about 10 percent. In single-mode case, the relief operation is not completed and unsatisfied demand still exists until the end of the planning horizon (minute 1440). However, by adding the second transportation mode, we are able to satisfy all demands by minute 1290 which is 2.5 hours before the end of the planning horizon. 145 Table 6.3 Optimization results for single and multimodal cases Case Description Objective Function Locations Selected Time of Last UD (min) CPU time (sec) Single Modal 38848500 4,5,8,10,12,13,17,19 1440 27 Multimodal 35110500 4,5,8,10,12,13,17,19 1290 3263 The most important comparison between the two cases is related to the CPU time to find the optimal solution in each case. Single modal formulation is solved to optimality (MIP gap = 0.1%) in only 27 seconds but it takes much longer (3263 seconds) to solve the multimodal numerical case. It is very interesting to notice that the number of variables and constraint in multimodal case is only about 10% more than single-mode case but it is much more difficult to solve the multimodal problem and it takes about 120 times longer to find the optimal solution in the second case. Figure 6.3 compares the performance of the relief operations for both cases over time. Total unsatisfied demand for both cases is shown for the duration of the operation. It can be seen that the two cases perform similarly for the first 8 hours of the operations; in fact no commodity is delivered to demand points during this time. However, the multimodal system has performed constantly better than the single-modal case after that initial period. 146 0 100000 200000 300000 400000 500000 600000 700000 0 240 480 720 960 1200 1440 To tal U ns ati sfi ed D em and operations time (min) Multi-Modal Single-Modal Figure 6.3 Comparison of performance for single and multi-modal cases It was explained that air transportation mode only covers the federal level network nodes. In other words, both modes cover all the nodes in federal level network but state level networks are covered only by the ground transportation mode. Analyzing the flow of commodities in the second case shows that for the federal network, 789500 units of commodities are sent by ground transportation mode versus 417000 units that are sent using the air mode. The market share of air transport is about 35 percent which shows the importance of fast transportation modes such as airplanes in emergency operations. It is important to notice that the number of available planes is one third of the number of trucks and the capacity of one plane is about 75% of the capacity of a trailer truck. 6.2.3 Analysis of Time-Step t Another parameter that affects the structure and behavior of the model is the length of time-step t. Time-step t is the length of time between two consecutive states that the problem is modeled. Selection of appropriate time-step 147 is a very important factor that can affect the performance of time-space networks dramatically. For each time period in the planning horizon, one layer of physical network will be added to the problem. This makes the problem size grow extremely fast with the number of time-steps in the planning horizon. For example if the planning horizon is 1 day, with the choice of time-step t = 5 minutes, 24 * 60 / 5 = 288 layers of the network is required to cover 1 day of operations. So to keep the problem at a reasonable size, it is favorable to have longer time-steps. On the other hand when t is short, the situation on the ground can be modeled in greater details which would not be possible with longer time-steps. For example if the time-step is 1 hour, it is only possible to model the state of the system at every hour and not at the times in between. So from the accuracy perspective, it is favorable to have shorter time-steps. Finding a reasonable time-step is an important modeling challenge. In selecting the time-step one should consider the level of accuracy that is required for that specific application and also the computational power that is available to them. In this section, it is tried to create and test numerical experiments with different values of time-step t and then analyze the findings in order to provide insight for other researchers or practitioners. In addition to computational aspects, changing the time-step length also affects some other elements of generating the mathematical formulation: 1. Link travel times 2. Capacity constraints 148 3. Objective function unit First, all link travel times for all transportation modes should be a multiplier of time-step t. For example when t = 5 minutes, travel time of a link cannot be 28 minutes or 47 minutes but it should be rounded to 30 and 45 minutes respectively . When the time-step is changed in a problem instance, for example from t = 5 to t = 10, a computer code is required to automatically recalculate the travel times and round them to the new time step. Second, loading and unloading capacity constraints at facilities are also required to be adjusted for different time-steps. For example, loading capacity of a given facility when t = 5 is equal to 1000 which means up to 1000 units of commodities can be loaded in that facility during that 5 minute time interval. If time-step lenght is changed to say t = 15 minutes, then it is necessary to adjust that loading capacity to 1000 * (15/5) = 3000 units. Third, the objective function in the proposed model is to minimize the sum of all unsatisfied demands and this summation is taken over all time periods. In the case that planning horizon is divided into N time periods of length t, the objective function summation involves N sets of variables. However, if a longer t is selected which is twice the previous t, then the same planning horizon consists of N/2 time periods, and the objective function value of these two cases will not have the same unit. To deal with this issue, it is recommended to normalize the unsatisfied demand over time and pay attention to the units that are used for those variables. 149 To analyze the effect of different values of t, 4 cases with values of t = 5, 10, 15, and 30 minutes are formulated and solved. Table 6.4 summarizes the problem size for these 4 cases. It is evident that the length of time step t has a huge impact on the number of variables and the number of constraint in the formulation. Table 6.4 Problem sizes for different values of time-step t Case Description Number of Variables Number of Constraints None-Zero Coefficients Input file size (Kb) t = 5 297194 90738 986268 19745 t = 10 157226 50125 525753 10456 t = 15 110570 36585 372297 7355 t = 30 63914 23033 218777 4257 Optimization results for these 4 cases are presented in Table 6.5. The stopping criteria for optimization code is set to be 0.5% optimality gap and objective function values are normalized to be of the same unit and comparable among all 4 cases. Table 6.5 Optimization results for different time-step length Case Description Objective Function Locations Selected Time of Last UD (min) CPU time (sec) t = 5 35155100 4,5,8,10,12,13,17,19 1295 888.39 t = 10 35238100 4,5,8,10,12,14,17,19 1290 124.34 t = 15 35199000 4,5,8,10,12,13,17,19 1290 79.16 t = 30 35344600 4,5,8,10,12,13,17,19 1290 62.45 150 Selecting the time-step length is a trade-off between modeling accuracy and solution time. Figure 6.4 illustrates the variations of the problem size and CPU time for the four cases modeled in this section. It is evident that when the length of the time-step is increased, the problem size and solution times both decrease. For this problem, t = 15 minutes seems to be a good trade-off between the accuracy and the problem size. For any specific application, it is recommended to initially perform similar analysis and then select the appropriate time-step length based on the required accuracy and availability of computational resources. 0 100 200 300 400 500 600 700 800 900 1000 0 50000 100000 150000 200000 250000 300000 350000 t = 5 t = 10 t = 15 t = 30 Nu mb er of Va ria ble s o r C on str ain ts Time Step Length Number of Variables Number of Constraints CPU time (sec) se ce co nd Figure 6.4 Variations of problem size and CPU time for different time steps 6.3 Sensitivity Analysis of the Main Input Parameters In this section, the main input parameters of the model are classified into 3 major categories: parameters of the facility location problem, parameters of the vehicle routing problem and parameters of the commodity flow problem. 151 Sensitivity analyses of the parameters in each of these categories are provided in the following subsections. 6.3.1 Sensitivity Analysis for Parameters of the Facility Location Number, location and capacity of the facilities in the network can have a major affect on the emergency response operations. From a list of potential locations, only a subset of them can be selected due to limitations of cost, equipment, or personnel. Also the capacities of each facility type can affect the response operations. In this section, the effect of variation of maximum number of facilities as well as loading capacity and unloading capacity at each facility is investigated. Maximum Number of Temporary Facilities As described earlier in the problem statement, there are three types of temporary facilities in FEMA?s supply chain structure. First, the mobilization centers (MOB) that receive the relief commodities from the permanent sources and forward them to federal operational staging areas (FOSA). Second, the FOSA that receive the commodities from permanent sources as well as mobilization centers and forward them to the state level facilities called state staging areas (SSA). Finally the state staging areas that receive flow from the FOSA and send them toward the final points of distribution (POD). In the current numerical example, 4 potential sites for MOBs are planned. Opening all facilities for operation provides the maximum capacity for relief operations, however; it might not always be possible to use all 4 facilities due to the high cost or limitations of the equipment and personnel. In order to investigate 152 the effect of this limitation, 5 numerical cases are formulated and solved. Table 6.6 introduces these cases and reports the optimization results. For these cases, it is assumed that all temporary facilities at lower levels (i.e. FOSAs and SSAs) are forced to be open. Table 6.6 Analysis of Number of Mobilization Centers Case No. Case Description Objective Function MOB Selected (Node number) 1 No MOB 3.78945 E+7 NA 2 1 MOB 3.55320 E+7 5 3 2 MOB 3.44932 E+7 4 , 5 4 3 MOB 3.44670 E+7 4 , 5 , 6 5 4 MOB 3.44580 E+7 4 , 5 , 6 , 7 In case 1, there is no MOB selected. It is shown that the relief operations can still proceed even without a mobilization center. The reason for that is the special structure of FEMA?s supply chain. As shown in figure 1.5, relief commodities can be send from the logistics centers and commercial storage sites directly to the federal operational staging areas and from there to each state and local area. However, having MOBs can provide more options and facilitate the flow of commodities and vehicles to the lower level facilities. Comparing case 1 and 2, by only opening one MOB in case 2, the objective function is considerably reduced. Figure 6.5 illustrates the effect of the number of MOBs on the objective function. 153 From the results presented in Table 6.6 and Figure 6.5, it can be concluded that adding more mobilization centers is beneficial in reducing the total unsatisfied demand. In Table 6.6, the objective function is constantly reduced when more MOBs become available. However, as shown in Figure 6.5, the improvements in objective functions become marginal when more that 2 MOBs are selected. Consequently, it is suggested to have a maximum of 2 mobilization centers in this specific numerical example. Figure 6.5 Effect of the number of mobilization centers on the objective function The same analysis is performed for the number of federal operational staging areas (FOSA). In this part, it is assumed that the 2 MOB are open and also all lower level facilities (SSA) are open. In order to investigate the effect of the number of FOSAs, 5 numerical cases are formulated and solved. Table 6.7 introduces these cases and reports the optimization results. In case 1 when there is no FOSA selected, the objective function is very high. In fact, no commodity is delivered in case 1 because without any FOSA the federal level and state level networks are disconnected. Based on the assumptions 3.4 3.45 3.5 3.55 3.6 3.65 3.7 3.75 3.8 3.85 0 1 2 3 4 Ob jec tiv e F unc tio n ( E+ 7) No. of MOB 154 of the FEMA?s supply chain structure, the commodities can only pass through FOSA nodes in order to be delivered to the state and local facilities. Figure 6.6 demonstrates the effect of the various number of FOSAs on the objective function. Having more FOSA has a positive effect in minimizing the objective function; however, this effect becomes very marginal for more than 2 FOSAs. Table 6.7 Analysis of the number of the federal operational staging areas Case No. Case Description Objective Function FOSA Selected (Node number) 1 No FOSA 5.82000 E+7 NA 2 1 FOSA 4.52235 E+7 8 3 2 FOSA 3.50445 E+7 8 , 10 4 3 FOSA 3.46125 E+7 8 , 9 , 10 5 4 FOSA 3.44932 E+7 8 , 9 , 10 , 11 2 2.5 3 3.5 4 4.5 5 5.5 6 0 1 2 3 4 Ob jec tiv e F un cti on (E +7) No. of FOSA Figure 6.6 Effect of the number of the federal operational staging areas on the objective function value 155 The next level of temporary facilities in the FEMA?s structure is the state staging areas (SSA). A similar study is conducted for the number of SSAs. There is a total of ten potential locations for SSAs in this problem. Number of MOBs and FOSAs are limited to 2 facilities for each type. Table 6.8 introduces these cases and summarizes the optimization results. Figure 6.7 illustrates the variation of objective function in each case. The objective function is constantly reduced in the first 5 cases but becomes steady after that. It can be concluded that for this numerical example, opening more than 4 state staging areas (2 in each disaster area) does not improve the performance of the operations. Table 6.8 Analysis of the number of the state staging areas Case No. Case Description Objective Function SSA Selected (Node number) 1 No SSA 5.82000 E+7 NA 2 max 1 SSA 4.66815 E+7 13 3 max 2 SSAs 3.81090 E+7 13,17 4 max 3 SSAs 3.56775 E+7 13,17,19 5 max 4 SSAs 3.50445 E+7 12,13,17,19 6 max 5 SSAs 3.50445 E+7 12,13,14,17,19 7 max 6 SAAs 3.50445 E+7 12,13,14,17,19,20 8 max 7 SSAs 3.50445 E+7 12,13,14,17,19,20,21 9 max 8 SSAs 3.50445 E+7 12,13,14,16,17,19,20,21 10 max 9 SSAs 3.50445 E+7 12,13,14,15,16,17,19,20,21 11 max 10 SSAs 3.50445 E+7 12,13,14,15,16,17,18,19,20,21 156 2 2.5 3 3.5 4 4.5 5 5.5 6 6.5 0 2 4 6 8 10 Ob jec tiv e F un cti on No. of SSA Figure 6.7 Effect of the number of state staging areas on the objective function value In the cases introduced in tables 6.6 to 6.8, the facility location constraints in each tier are treated independently. In fact, the problem is first solved to find the optimal number of MOBs then it is solved for FOSAs and then SSAs respectively. However, when the resources (e.g. equipment and personnel) can be shared among the different facility types, the maximum number of each facility type is not independent anymore. For example instead of having 2 MOBs and 2 FOSAs, it might be beneficial to have 1 MOB and 3 FOSAs. Nineteen cases are generated and solved for the maximum number of facilities of all types. Table 6.9 lists these cases and optimization results. In case 0, no temporary facility is selected and no commodity can be delivered because the supply chain is disconnected. The relative objective function column in table 6.9 shows the ratio of the objective function of each case to the maximum objective function value given in case 0. 157 Figure 6.8 illustrates the objective function values and optimization CPU times for the cases introduced in table 6.9. Table 6.9 Sensitivity analysis of the total number of the temporary facilities Case Max No of Facilities Objective Function Relative O.F. (%) Locations Selected CPU Time(sec) 0 0 5.82 100 - 0.3 1 1 5.82 100 4 37 2 2 4.66815 80.21 8,13 137 3 3 4.55725 78.30 8,12,13 423 4 4 4.1403 71.14 8,10,13,17 650 5 5 3.89115 66.86 4,5,10,13,17 364 6 6 3.70395 63.64 5,8,10,13,17,19 478 7 7 3.56775 61.30 4,5,8,10,13,17,19 257 8 8 3.50445 60.21 4,5,8,10,12,13,17,19 133 9 9 3.48225 59.83 4,5,8,10,11,12,13,17,19 190 10 10 3.46305 59.50 4,5,8,10,11,12,13,17,19,20 115 11 11 3.4566 59.39 4,5,8,9,10,11,12,13,17,19, 20 128 12 12 3.450825 59.29 4,5,8,9,10,11,12,13,15,17, 19,20 156 13 13 3.4485 59.25 4,5,6,8,9,10,11,12,13,15,17,19,20 52 14 14 3.447 59.23 4,5,6,8,9,10,11,12,13,14,15,17,19,20 60 15 15 3.4461 59.21 4,5,6,7,8,9,10,11,12,13,14, 15,17,19,20 24 16 16 3.4458 59.20 4,5,6,7,8,9,10,11,12,13,14, 15,17,19,20,21 6 17 17 3.4458 59.20 4,5,6,7,8,9,10,11,12,13,14, 15,17,18,19,20,21 5 18 18 3.4458 59.20 4,5,6,7,8,9,10,11,12,13,14, 15,16,17,18,19,20,21 6 158 0 100 200 300 400 500 600 700 40 50 60 70 80 90 100 0 2 4 6 8 10 12 14 16 18 Un sat isf ied De ma nd (%) No. of Temporary Facilities Objective Function CPU-Time CP U T ime (s ec) Figure 6.8 Objective function and CPU time versus the number of temporary facilities Analyzing the results shows that increasing the number of the available temporary facilities can effectively help the operations and reduce the total unsatisfied demand. However, after a certain number of facilities, adding more facilities does not reduce the objective function value. For example in current problem, selecting more than 8 temporary facilities has very marginal benefits. The CPU time is relatively low for the first 3 cases because the combinations are limited. CPU times considerably increase for case 3 to 7 mainly due to the increase in the number of potential combinations. For the final 6 cases, CPU time is reduced again because ample capacity is available and opening more facilities does not affect the objective function and does not change the flow of commodities or vehicles. Loading and Unloading Capacities Any facility that sends or receives the relief commodities is subject to a limited capacity for loading and unloading (for mathematical formulations refer to 159 section 3.5.5). These constraints are mainly due to the limitations of equipment or personnel during the emergency response operations. In this section, it is intended to evaluate the effect of these capacities on the performance of the model and the optimization results. Table 6.10 lists loading and unloading capacities (cft/hr) for each transportation mode and each facility type that are used in the base case scenario. Table 6.10 The loading and unloading capacities Ground Transportation Air Transportation Facility Type Loading Unloading Loading Unloading LC - CSS - VEN 60000 60000 18000 18000 MOB 48000 48000 9000 9000 FOSA 36000 36000 - - SSA 24000 24000 - - POD 12000 12000 - - To evaluate the effect of variations in loading capacity, five different values are considered each with 50%, 75%, 100%, 125% and 150% of the original value. The same five variations are also considered for unloading capacity. Also to evaluate the joint effect of loading and unloading capacity on the performance of the model a matrix is generated to consider all 25 combinations. Table 6.11 reports the optimization results for these 25 cases. 160 Table 6.11 Optimization results for various Loading and Unloading Capacities Objective Function Unloading Capacity 50% 75% 100% 125% 150% Loading Capacity 50% 4.19625 4.19235 4.19235 4.19235 4.19235 75% 4.19235 3.7632 3.69975 3.66765 3.66765 100% 4.19235 3.7578 3.50445 3.468675 3.44715 125% 4.19235 3.7578 3.49965 3.3771 3.35835 150% 4.19235 3.7578 3.49965 3.37185 3.31095 Focusing on the 3rd row it is shown that when loading capacity is kept at 100%, increase in unloading capacity from 50% to 150% constantly improved the objective function value. The same behavior is shown in the 3rd column for original unloading capacity when loading capacity varies from 50% to 150%. However, comparing the objective function values of the first row, it is evident that when loading capacity is fixed at 50%, the problem is not sensitive to the unloading capacity anymore. The similar behavior is observed in the first column when unloading capacity is fixed at 50% and extra loading capacity has no benefits. It is concluded that extra capacity at facilities can be useful in reducing the objective function. However, these additional capacities are beneficial only when both loading and unloading capacities are increased at the same time and proportionally. If one capacity is kept considerably low, additional capacity of the other type is not effective anymore. To visualize, figure 6.9 illustrates the 161 objective function values versus different values of the loading and unloading capacities. Ob ject ive Fun ctio n V alue 4.4-4.5 4.2-4.4 4-4.2 3.8-4 3.6-3.8 3.4-3.6 3.2-3.4 3-3.2 Figure 6.9 Objective function value versus variations in the loading and unloading capacity 6.3.2 Sensitivity analysis for Parameters of Vehicle Routing Vehicle routing problem is a major part of the proposed integrated logistics model. Variations in the inputs of vehicle routing can drastically impact the behavior and results of the entire model. In this section, a series of analysis is performed to investigate the nature and extent of these effects. In the following, sensitivity analysis is performed on the number of vehicles of each type, capacity of vehicles, and travel speeds of the vehicles. Number of Available Vehicles Having more vehicles is always favorable from operator?s perspective because it can provide more capacity for faster and easier delivery of relief items. However, the number of available vehicles is limited especially during the initial 162 periods of disaster response operations. Analysis of the effects of different number of vehicles can provide invaluable insight to the problem for planning purposes. In table 6.12 ten cases are tested with various numbers of vehicles for each mode. For all these cases, the vehicles of ground transportation mode are 53 ft tractor trailers with 6000 cft capacity. At the beginning of the operations, trucks are evenly distributed among 3 source nodes (LG, CSS, VEN). The vehicles of air transportation mode are C-130 Hercules cargo planes with 4500 cft capacity that are located at CSS and VEN nodes at the beginning of the operations. Table 6.12 Sensitivity analysis of number of vehicles Case No. No. of Trucks No. of Planes Objective Function (E+7) CPU-Time (sec) 1 12 4 4.1345 68.9 2 24 8 4.08015 372.42 3 36 12 3.82065 237 4 48 16 3.63045 848 5 60 20 3.51075 610 6 72 24 3.44025 50 7 84 28 3.39105 43.45 8 96 32 3.3546 4.86 9 108 36 3.32715 8.05 10 120 40 3.32565 8.72 In cases 1 to 10, the numbers of vehicles of both modes are gradually increased. Case 5 is similar to the base case in previous subsections. CPLEX commercial solver is used to optimize these case studies. Stopping criteria is set to be 0.1% optimality gap. Figure 6.10 illustrates the objective function and the CPU times for optimal solutions in each case. 163 0 100 200 300 400 500 600 700 800 900 3 3.2 3.4 3.6 3.8 4 4.2 1 2 3 4 5 6 7 8 9 10 Ob jec tiv e F un cti on Case Number Objective Function CPU-Time CPU -T im e Figure 6.10 Effect of the number of vehicles on the objective function and the CPU time The value of objective function decreases as the number of vehicles increase. It complies with expectations since more vehicles provide higher capacity and faster delivery of relief items which minimizes the objective function value. The rate of decrease in objective function is faster for the first five cases compared to the last five cases. For example, the objective function value decreases by 17% from case 1 to case 5. The same measure is only about 6% from case 5 to case 10. It can be concluded that when number of vehicles are increased in cases 6 to 10, the other constraints become binding. In that case, it is recommended to invest in other parts of the system and increase other capacities such as loading and unloading capacities at transfer facilities. CPU-time variations are very interesting. At the beginning, the CPU time is relatively low because there is a very small number of vehicles available and there is not much room for improvement. As the number of vehicles increase, the combinations for vehicle routing problem increases rapidly. As a result, the CPU time to find the optimal solution is grown considerably. At the end, it might be 164 surprising to notice that the CPU time is decreased. The reason is that for cases 6 to 10, a large number of vehicles are available. When there are ample vehicles available, the model can easily assign a new vehicle from depot to any required task. On the other hand, when the number of vehicles are limited, a great deal of time is spend to search all possible combinations and make sure to assign the vehicles to the best possible task. It is concluded that the number of vehicles not only affects the model?s results such as objective function value, but it also affects the difficulty of solving the model and CPU time to find the optimal solution. When the number of vehicles is very low or very high, it is much easier and faster to solve the model. For an in-between range of vehicle numbers, it can become very difficult and time consuming to find the optimal solution. This range is problem-specific and can depend on the other model inputs as well. Researchers and practitioners should be aware of this behavior and perform the similar analysis for a range of vehicle numbers that is specific to their specific application. Capacity of the Vehicles Another factor that can affect the performance of the entire model is the capacity or type of vehicles that are used in the response operations. The general conception is that higher capacity is always better. That might be true; however it might not be possible to always use the largest vehicle in the fleet. In this section it is intended to analyze the effects of having vehicles with different capacities. Ten different cases are tested. In the first 5 cases, the capacity of ground transportation vehicles are changed from 2000 cft to 10000 cft while the capacity 165 of planes are fixed at 4500 cft. Then in cases 6 through 10, the capacity of air transportation vehicles are changed from 1500 cft to 7500 cft while capacity of truck are kept fixed at 6000 cft. Table 6.13 lists these cases and reports the optimization results. Table 6.13 Sensitivity analysis of the vehicle's capacity Case No. Capacity of Trucks Capacity of Planes Objective Function (E+7) CPU-Time (sec) 1 2000 4500 4.53145 13 2 4000 4500 3.86310 1333 3 6000 4500 3.51075 610 4 8000 4500 3.42365 532 5 10000 4500 3.35250 285 6 6000 1500 3.62205 258 7 6000 3000 3.55590 148 8 6000 4500 3.51075 610 9 6000 6000 3.50175 647 10 6000 7500 3.51075 665 Figure 6.11 illustrates the variations of the objective function value and the CPU-time for these 10 cases. For cases 1 to 5 and 6 to 10, the objective function is decreased when the vehicle capacity is increased. Comparing the first 5 cases (left side of the graph) to the last 5 cases (right side) it is evident that the slope of the objective function is steeper for the left curve. It means that the problem is more sensitive to the capacity of trucks rather than to the capacity of planes. This might be simply due to the fact that the ground transportation does the majority of the deliveries in this problem. For example in the base case (case 3 166 and 8), total amount of commodity-miles that are shipped with ground transportation is about 250 million cft.mile versus about 107 million cft.mile for air transportation. Share of ground transportation in this problem is about 70%. Figure 6.11 Sensitivity analysis of the vehicle's capacity Travel speed The other factor in the vehicle routing problem that can affect the response operations is the travel time between the nodes. Faster vehicles are favorable from two perspectives. First, the flow of relief commodities through the network can happen faster. Second, empty vehicles can travel faster and reach the pick-up nodes to start another round of deliveries in a shorter period of time. On the other hand, in some cases the travel speeds might be reduced due to the changes in disaster environment such as inclement weather or flooding. In order to investigate the effects of travel speed on the operation?s performance, 8 different cases are generated and tested. Table 6.14 lists these cases and summarizes the optimization results. 0 200 400 600 800 1000 1200 1400 3 3.2 3.4 3.6 3.8 4 4.2 4.4 4.6 4.8 1 2 3 4 5 6 7 8 9 10 Ob ject ive Func tion Case Number Objective Function Objective Function CPU-Time CPU-Time CPU -Ti me (s ec) 167 The first case is the base case scenario. In the base case, it is assumed that average travel speed for ground transportation is 50 mph for the arcs in the federal level network and 40 mph for the arcs of the state level network. Average air travel speed is also assumed to be 200 mph. Table 6.14 Sensitivity Analysis for network Travel Speeds Ground Travel Speed Case No. Federal Network State Network Air Travel Speed Objective Function (E+7) CPU-Time (sec) 1 50 40 200 3.51075 610 2 60 48 240 3.06945 54.19 3 40 32 160 4.18710 3116 4 50 40 250 3.48735 678 5 50 40 300 3.47805 3519 6 50 40 400 3.46155 1024 7 40 40 200 4.06125 357 8 50 30 200 3.6528 961 In case 2 and case 3, it is assumed that for both transportation modes and the entire network the travel speed is increased by 20% in case 2 and decreased by 20% in case 3. It can be seen that if travel speed is increased by 20%, the objective function can be improved by 12.6%. However, if the travel speed is decreased by 20%, then the objective function is increased by 19.3%. It is evident that travel speed has a major effect on the efficiency of the operations. Comparing the first 3 cases, it can be concluded that faster transportation improves the performance but having slower transportation can have a larger negative impact. Case 4, case 5 and case 6 are similar to base case but the average air travel speed is increased from 200 mph to 250, 300 and 400 mph respectively (25%, 168 50%, and 100% increase). In these cases, the objective function is improved but only by 0.7 percent, 0.9 percent and 1.4 percent, respectively. It can be concluded that using faster planes can improve the operations but this improvement is very marginal when the ground transportation is not enhanced proportionally. In case 7, it is assumed that the ground travel speed in federal level network is reduced from 50 mph to 40 mph while other inputs are similar to the base case. In case 8, only the speed in state level networks is reduced from 40 mph to 30 mph and all other inputs are unchanged. In case 7, objective function is higher than base case by about 16% but the same measure is only 4% higher for case 8. Comparing case 7 and 8, it seems that travel speed at federal level network has a stronger effect than travel speed at the state level networks. In other words, if we can choose to improve the conditions of the roads in either federal level or state level roadways, it can be more rewarding to improve the federal level links. 6.3.3 Sensitivity analysis for Parameters of the Commodity Flow In this section, sensitivity analysis is performed to investigate the effect of various amounts of supply, random demands and the relative urgency factor. Sensitivity to Supply In the base case scenario, the entire supply for one day of operations is assumed to be available at the beginning of the operations. Total supply for one day of operations is 600,000 units that include 400,000 units of commodity 1 and 200,000 units of commodity 2. At the beginning of the operations, supply is stored at three main source facilities. 40% of supply is stored at the logistics center site in Atlanta, GA. 20% is stored at commercial storage site in Charlotte, 169 NC and 40% of supply is stored at the vendors site in Nashville, TN. The availability of supply can play a major role in disaster response operations. In this section, 8 different cases are generated and tested to evaluate the effect of availability of different levels of supply. Table 6.15 introduces these cases and summarizes the optimization results. Table 6.15 Sensitivity analysis for supply availability Case No. Supply Availability Objective Function (E+7) CPU-Time (sec) Locations Used 1 50% 4.23225 11 4,5,8,10,12,13,17,19 2 60% 4.02015 8.75 4,5,8,10,12,13,17,19 3 70% 3.84705 5.7 4,5,8,10,12,13,17,19 4 80% 3.70785 4.39 4,5,8,10,12,13,17,19 5 90% 3.59085 5.84 4,5,8,10,12,13,17,19 6 100% 3.50445 7.78 4,5,8,10,12,13,17,19 7 110% 3.49905 7.08 4,5,8,10,12,13,17,19 8 120% 3.49905 8.97 4,5,8,10,12,13,17,19 Figure 6.12 illustrates the results of table 6.15. Case 6 is the base case in which supply is equal to 100% of the demand and can be used as a benchmark. From case 5 to case 1, the amounts of supplies available at the sources are gradually reduced. It can be seen that shortage of supplies can strongly affect the results. For example in case 1, shortage of supplies has resulted in about 20% higher objective function value. On the other hand, in cases 7 and 8 there are extra supplies available at each source node. It is evident that having additional supplies has a marginal effect on improving the objective function (less than 1% reduction). The reason is that in these cases the other constraints of the problem 170 become binding. In fact, limitations such as the transportation capacity and the facility constraints limit the amounts of supply that can be delivered and having extra supply cannot help to reduce the objective function. 0 5 10 15 20 25 30 3 3.2 3.4 3.6 3.8 4 4.2 4.4 1 2 3 4 5 6 7 8 Ob jec tiv e F un cti on Case Number Objective Function CPU-Time CPU -T im e ( sec) Figure 6.12 Sensitivity analysis for amount of supplies available Sensitivity to Demand Locations and amounts of demands are other variables that affect the details of the response operation. In this section, 10 cases with random demand values are generated and solved to optimality. The demands in current numerical example are located at 20 points of distribution (POD) that are spread over 2 disaster areas. Random populations are assigned to each POD and demand for each commodity is generated based on the population of each POD. The total population of 2 disaster areas are kept fixed but population of PODs are different in each of the 10 cases. Random population for each POD is generated using a uniform distribution between 300 and 1700 people. 171 Table 6.16 and figure 6.13 represent the optimization results. Even though the demands at each POD are different, the values of the objective function across the cases are not that different. The coefficient of variation for objective function value among the 10 random cases is relatively small (only 2.3%) which shows that the proposed model is successful in managing the demand variations. Also, it is important to notice that the same set of temporary facility locations are selected for all ten cases. This observation is favorable and can be interpreted as a good measure of robustness of the model in case of fluctuations in demand. Table 6.16 Sensitivity analysis for variations in demand Case No. Objective Function (E+7) CPU-Time (sec) Locations Used 1 3.50445 7.78 4,5,8,10,12,13,17,19 2 3.46545 12.84 4,5,8,10,12,13,17,19 3 3.5199 11.31 4,5,8,10,12,13,17,19 4 3.53715 10.55 4,5,8,10,12,13,17,19 5 3.49095 16.45 4,5,8,10,12,13,17,19 6 3.48825 7.7 4,5,8,10,12,13,17,19 7 3.73995 19.91 4,5,8,10,12,13,17,19 8 3.4392 6.09 4,5,8,10,12,13,17,19 9 3.50805 6.55 4,5,8,10,12,13,17,19 10 3.5262 18.66 4,5,8,10,12,13,17,19 Average 3.521955 11.784 Coef.of var. 0.0232 0.429 172 0 30 60 2 2.5 3 3.5 4 4.5 5 1 2 3 4 5 6 7 8 9 10 Ob jec tiv e F un cti on (E +7) Case Number Objective Function CPU-Time CPU -T im e ( sec) Figure 6.13 Optimization results for random demand values Relative Urgency factor As defined in section 3.5.2, citRU is relative urgency of one unit of commodity c, in node i at time t. It is a weight factor in the objective function to enforce the importance of one commodity over another or one demand point over another demand point. For example, if one unit of commodity 1 is more important than 1 unit of commodity 2, then RU for commodity 1 should be higher than RU for commodity 2. In this section, effect of different values of the relative urgency factor is investigated. First, the effect of using different relative urgency factors among commodities is investigated. There are 2 commodities in the current numerical example. Five cases are generated to test the different combinations of weight for these 2 commodities. In the base case, RU is equal to 1 for all commodities. In cases 2 and 3, the priority is given to the 1st commodity (C1). In cases 4 and 5, the 173 priority is given to the 2nd commodity (C2). Table 6.17 introduces these cases and reports the optimization results. In Table 6.17, UR1it and UR2it are the relative urgency factors for commodity 1 and 2. The two columns on the right side of the table represent the sums of unsatisfied demand of each commodity by itself. Total unsatisfied demand is the sum of these two columns. Table 6.17 Sensitivity analysis of the relative urgency factor for each commodity Case No. UR1it UR2it Total Unsatisfied Demand 1 1 1 35044500 22891500 12153000 2 2 1 35044500 20714000 14330500 3 5 1 35044500 20714000 14330500 4 1 2 35069500 26216000 8853500 5 1 5 35096500 26264000 8832500 Figure 6.14 better illustrates the results of Table 6.17. In cases 2 and 3, urgency factor of commodity 1 is increased to 2 and 5. As a result, total sum of unsatisfied demands of commodity 1 is decreased compared to the base case (case 1). On the other hand, the same measure for commodity 2 is increased in cases 2 and 3. It can be said that, while the total objective function is the same, the shipments of commodity 1 are delivered faster (in earlier times) than commodity 2 because of the higher urgency factor. In cases 4 and 5, the relative urgency factor for commodity 2 is increased. In these cases, demands for commodity 2 are satisfied earlier and share of 174 unsatisfied demand of commodity 2 in the objective function is decreased compared to the base case scenario. 22.8915 20.714 20.714 26.216 26.264 12.153 14.3305 14.3305 8.8535 8.8325 0 5 10 15 20 25 30 35 40 1 2 3 4 5 Mi llio ns SUM C1 SUM C2 Un sat isf ied De ma nd Case Number Figure 6.14 Sensitivity analysis of the relative urgency factor for each commodity Figure 6.15 can help to better understand the effect of the RU factor on the performance of the response operations. In Figure 6.15, variations of unsatisfied demand for each commodity over time are shown for cases 1, 2 and 4. The solid line shows commodity 1 and the dashed line belongs to commodity 2. In case 1 (left-side graph), at time zero the demand for commodity one is 400,000 units versus 200,000 units for commodity 2. RU is equal to 1 for both commodities. It can be seen that over time, the demands for both commodities are satisfied almost proportionally and by the end of the operations both lines get to zero almost at the same time. In case 2 (Figure 6.15, center graph), the priority is given to commodity 1. Even though the initial unsatisfied demand for commodity 1 is much higher, as 175 the time goes by, the demands of commodity 1 are satisfied more rapidly. In fact, for the first 13 hours of the operation all of the deliveries are the shipments of commodity 1 only. The dashed line shows the rate of demand satisfaction for commodity 2. Compared to the left side graph, the demands of commodity 2 are satisfied much later in Case 2 due to the higher priority of commodity 1. On the contrary in case 4 (Figure 6.15, right-side graph), the priority is given to commodity 2. Consequently, the demands of commodity 2 (dashed line) are satisfied much earlier compared to case 1 or case 2. Figure 6.15 Effect of relative urgency factor on variations of unsatisfied demand over time Overall, it is shown that using relative urgency factors can be very effective when there is a reason to give priority to one commodity over another. Usually, high priority commodities are needed in much smaller quantities compared to the commodities with normal priority. For example demand for medical supplies are usually in much lower volumes compared to clothing or construction items but with a much higher priority. When there is no urgency factor in the model, the commodities with higher demand volumes tend to be given the priority in order to minimize their demand. However, when there is a small amount of demand for a commodity with high priority (e.g. medical 0 50 100 150 200 250 300 350 400 450 0 240 480 720 960 12001440 To tal Uns ati sif ied d em and operations time (min) Case1-C1 Case1-C2 0 50 100 150 200 250 300 350 400 450 0 240 480 720 960 12001440 To tal U ns ati sif ied d em and operations time (min) Case2-C1 Case2-C2 0 50 100 150 200 250 300 350 400 450 0 240 480 720 960 12001440 To tal U ns ati sif ied d em and operations time (min) Case4-C1 Case4-C2 176 supplies), it is very effective to use a higher RU factor in the model to make sure that the priority is given to that commodity regardless of its small quantity. Relative urgency factors can also be used to enforce priorities among different points of distribution. If the demands for relief commodities at one or more locations have priorities, higher RU factors for those locations should be used to imply these priorities. The effect of using higher RU factors to prioritize a subset of locations is investigated in the following. Using the same numerical experiment from previous section, 2 numerical cases are simulated and compared. Case 1 is the base case where ticRUcit ,,1 ?= . In Case 2, one demand node in each disaster area is selected to have a higher priority. In the state of Mississippi, the POD located at node 29 is selected and from demand nodes in the state of Louisiana, the POD located at node 41 is selected to have a higher priority. These 2 nodes are selected among all PODs because they were distant locations and in the results of the base case scenario, they were the 2 nodes that received relief items later than other nodes. For these 2 nodes, a relative urgency factor of 2 is applied to the demands of all commodities and all time periods. (For geographical locations of the nodes of the network refer to Table 4.1) Case 1 is simulated with ticRUcit ,,1 ?= and Case 2 is formulated and solved with tcRU c t ,129 ?= and tcRU c t ,141 ?= , and all other RU = 1 for the rest of the demand nodes. Figures 6.16 and 6.17 illustrate the optimization results of the both cases for POD nodes 29 and 41. In both figures, it is evident that the required priorities are successfully enforced by using higher RU factors. In figure 6.16, before enforcing priorities, the last demand is satisfied by t = 960 minute; 177 however, after using priorities the last demand is satisfied by t = 510 minute which is much faster. Figure 6.16 Effect of using higher priority for POD 29 In figure 6.17, the effect of using and not using priorities for POD node 41 is illustrated. In case 1, when no priority is required, the last demand is satisfied by time t = 1290 minutes. However, by using priority in Case 2, the time of the last unsatisfied demand is reduced to only t = 870 minutes. These analyses show that the relative urgency factors can be successfully used to give priority to the demands of certain PODs if required by the user. However, the user should be aware of the fact that assigning these priorities would only improve the demand satisfaction rates for the intended nodes. Meanwhile, since the limited resources are directed to high priority nodes, the demand nodes with lower priorities would perform worse than before. 0 2000 4000 6000 8000 10000 12000 0 240 480 720 960 1200 1440 Un sat isf ied De ma nd of PO D 29 Operations time (min) Case 1 Case 2 178 Figure 6.17 Effect of using higher priority for POD 41 6.4 Summary and Conclusions In this chapter, in-depth analyses of different aspects of the proposed mathematical model were presented. These analyses were provided in order to better illustrate the capabilities of the model and also examine model?s sensitivity in various circumstances. The analyses in this chapter were divided into two main categories: sensitivity analysis of the structural parameters of the model, and sensitivity analysis of the main input values of the model. It was shown that changing some input parameters not only affects the optimization results but it can also largely change the problem size and solution computation times. The structural parameters investigated in section 6.1 include the number of commodities C, the number of transportation modes M, and the length of time- step resolution t. To test the effect of the number of commodities, four cases are considered each with one, two, three and four commodities respectively. It is 0 5000 10000 15000 20000 25000 30000 35000 40000 45000 50000 0 240 480 720 960 1200 1440 Un sat isf ied De ma nd of PO D 41 Operations time (min) Case 1 Case 2 179 shown that the number of variables and the number of constraints increase linearly with the number of commodities. Number of variables grows faster than number of constraints. More importantly, it is shown that when number of commodities increase, the problem becomes much harder to solve and the computation times rapidly increase. This increase in computation time was exponential in the range tested. Number of the transportation modes is also an important factor. Using multimodal compared to single modal transportation not only increases the size of the formulation but also makes the problem much harder to solve. In multimodal systems there is the question of which transportation mode to utilize and how to balance the commodity flows among different modes. Another concern in multimodal environments is considering the intermodal transfers. From application?s perspective it is important to provide suitable facilities and required equipments in order to transfer the relief commodities between transportation modes quickly and efficiently. From modeling perspective it is important to consider the properties of each transportation mode and correctly model the delays during intermodal transfers. Multimodal problem is more difficult to solve. In fact, for the current numerical example the single modal formulation was solved in only 27 seconds but it took much longer (3263 seconds) to solve the multimodal numerical case. It is very interesting to notice that the number of variables and constraint in multimodal case was only about 10% more than single-mode case but it was 180 much more difficult to solve the multimodal problem and it took about 120 times longer to find the optimal solution in multimodal problem. The other structural parameter in the model is the length of time-steps t. Time-step t is the length of time between two consecutive states that the problem is modeled. Selection of appropriate time step is a very important factor that can affect the performance and accuracy of the time-space networks dramatically. When t is short, the situation on the ground can be modeled in greater details which would not be possible with using longer time-steps. So from accuracy perspective, it is favorable to have shorter time-steps. On the other hand, for each time period in the planning horizon, one layer of physical network will be added to the time-space structure of the problem. This makes the problem size grow extremely fast with the number of time-steps in the planning horizon. To keep the problem size manageable, it is preferred to have longer time-steps. Finding a reasonable time-step t is an important modeling challenge. In selecting the time-step, one should consider the level of accuracy that is required for that specific application and also the computational power that is available to them. For this problem, t = 15 minutes seemed to be a good trade-off between the accuracy and the problem size. For any specific application, it is recommended to initially perform similar analysis and then select the appropriate time-step length based on the required accuracy and availability of computational resources. The main input parameters investigated in section 6.2 are divided into 3 categories: parameters of the facility location problem, parameters of the vehicle routing problem and parameters of the commodity flow problem. 181 It is shown that the number, location and capacity of the facilities in the network can have a major affect on the emergency response operations. In analysis of the number of facilities, it was concluded that adding more facilities is beneficial in reducing the unsatisfied demand over time. However, these improvements became marginal when more that 2 MOBs, 2 FOSAs and 4SSAs were selected. The loading and unloading capacity in each facility was shown to impact the flow of commodities and how the low capacities would increase unsatisfied demand. Numerical studies tested a range of capacities for both loading and unloading capacity factors. It was concluded that investments to expand the capacity should improve both capacities at the same time. If one of the capacities is kept at a fixed level, then additional capacity for the other type remains unused. For vehicle routing problem, the number of available vehicles was one of the main factors. It was concluded that the number of vehicles not only affected the model?s results such as objective function value, but it also affected the difficulty of model and CPU time to find the optimal solution. When the number of vehicles was very low or very high, it was much easier and faster to solve the model. For an in-between range of vehicle numbers, it became very difficult and time consuming to find the optimal solution. This range is problem-specific and can depend on the other model inputs as well. Researchers and practitioners should be aware of this behavior and perform the similar analysis for a range of vehicle numbers that is appropriate for their specific application. 182 Another factor that affected the performance of the entire model was the capacity or type of vehicles used in the response operations. The general conception is that higher capacity is always better. However it might not always be possible to use the largest vehicle in the fleet. Analysis of the capacity of trucks and planes in the system indicated that the problem was more sensitive to the capacity of trucks versus the capacity of planes. This might be simply due to the fact that the ground transportation does the majority of the deliveries in this problem. The other factor in vehicle routing problem that affected the response operations was the travel speeds. Faster vehicles are favorable from two perspectives. First, the flow of relief commodities through the network can happen faster. Second, empty vehicles can travel faster and reach the pickup nodes to start another round of deliveries in a shorter period of time. However, the travel speed might be reduced in the disaster area due to the inclement weather or road blockings. It was shown that lower travel speeds would increase the unsatisfied demands over time. Comparing the travel speeds, it was shown that travel speeds at federal level network had a stronger effect than travel speed at the state level networks. In other words, if we could choose to improve the conditions of the roads in either federal level or state level roadways, it would be more rewarding to improve the conditions of the federal level links. In commodity flow problem, sensitivity analysis was performed to investigate the effect of various amounts of supply, random demands and the 183 relative urgency factor. In analyzing available supplies, it was shown that the lack or delay of supplies at source nodes had a large negative effect on the objective function value. On the other hand, when supplies were abundant, the objective function was improved only slightly. In fact, other limitations such as transportation capacity and facility constraints limited the amounts of supply that could be delivered and having extra supply could not help. Locations and amounts of demands are other factors that affect the details of the response operation. Variability in demand locations and amounts is a negative factor for emergency response operations. Ten cases with random demand values were generated and solved to optimality. It was shown that in spite of variations in demand, the proposed model was successful in finding good solutions. Equally important, the same set of facility locations were used for all random cases that could be interpreted as a good measure of robustness of the model in case of fluctuations in demand. Finally, the effect of relative urgency factor was tested. It was proven that using relative urgency factor can be very effective when there was a reason to give priority to one commodity over another or one demand node over another. The demands for commodities or nodes with higher RU factor were satisfied at very earlier times compared to the other commodities or nodes. Use of the RU factor is highly recommended when there is a small amount of demand for a commodity with high priority (e.g. medical supplies).It is very effective to use a higher RU factor in the model to make sure that the priority is given to that commodity/node regardless of its small quantity. 184 Chapter 7: Prepositioning and Equity In this chapter 2 new subjects are investigated that are very important in emergency response operations. In section 7.1, first the concept of prepositioning is introduced and then the mathematical formulation to model the prepositioning problem is described. After that, a set of numerical experiments are conducted to illustrate the potential benefits of considering prepositioning in emergency response operations. In section 7.2, equity constraints that were previously introduced in chapter 3 are further investigated then a new set of equations are proposed that can model the equity constraints more effectively. Numerical experiments illustrate the successful use of the new equity constraints. 7.1 Prepositioning of supplies and vehicles Planning and preparedness play a vital role in disaster management and emergency response. After a disaster strikes, the initial unavailability of supplies or the slow pace in mobilizing them can cause major delays in emergency response that would result in increased loss of life and human sufferings. Prepositioning is a valuable tool for emergency response organizations to enhance their emergency response capacity and preparedness for responding to large-scale disasters. However, effective prepositioning of supplies and personnel is not an easy task. Uncertainty about future disasters and also the high costs of inventory and maintenance are some of the obstacles in effective prepositioning for large- scale emergency operations. 185 A number of researches in recent years have considered the prepositioning of supplies for emergency operations (Rawls and Turnquist, 2010). These studies emphasized the importance and benefits of considering strategic inventory and prepositioning of supplies. In this section, it is tired to use the mathematical model introduced in chapter 3, in order to model and optimize the prepositioning of supplies as well as transportation capacity for disaster response. It should be emphasized that technically, prepositioning problem is usually considered at the planning or strategic level. Conversely, the mathematical model proposed in chapter 3 is a tool to model and optimize emergency response operations at the operational level. Nevertheless, the unique capabilities of the proposed model can still be used to solve the prepositioning problem. To do so, it is required to generate a wide range of potential disaster scenarios in advance and solve them with the proposed model and then implement the aggregated outcomes of all scenarios in planning or strategic level decision making. 7.1.1 Mathematical Formulation In the mathematical model proposed in chapter 3, the amount of exogenous supply for each commodity at each node of the network was assumed to be a parameter that was given for each problem instance. Equation (7.1) below is similar to Eq. (3.2) from chapter 3: tmciSupCXXT XCXXTX c it c ti m mmc kti j cm ttji c it m mcm it j cm ijt mm jim ,,,)1()( )( ?=?? ?++ ? ? ? ? ? ? ? ? ??? ? (7.1) 186 For the initial time period, ? ? j cm ttji jimX )( is equal to 0 as well as ? ? ? ? ? m mmc kti mmXT )( and c tiCX )1( ? . So, if the parameter citSup is being considered as a variable of the model for the initial time period (t = 0), then we have: 0,,0 =?=?+?? tciSupCXX citcit m j cm ijt (7.2) Equation 7.2 is the new constraint for conservation of flow for the initial time period. Equation 7.2 requires that for any given node and any given commodity, the sum of all commodities that are shipped from i to any node by any transportation mode plus the amount being stored at that node for the next time period minus the exogenous supply is equal to zero. Since the total supply available at time 0 is limited, we need to add a new constraint to the problem to limit the amount of total initial supply cSup0 : 0,0 =?=? tcSupSup c i c it (7.3) Based on FEMA?s recommendations, the supplies can also be prepositioned at the temporary facilities. However, a temporary facility must be open and have enough storage capacity to hold prepositioned supplies. Equation 7.4 enforces these requirements: 0, =????? tWiLocScapSup tiit c c it (7.4) By adding the equations 7.2, 7.3 and 7.4 to the main model, we can solve for the optimal amounts and location for prepositioning of the available supplies. Prepositioning is not limited to the relief commodities but can be considered for the vehicles that transport these commodities. Finding the optimal 187 location for vehicles of each transportation mode for the initial time period is also very important. Using the concepts similar to the ones used for prepositioning of supplies, we will have: 0,,0 =??=?+? tmNiAVCYY mitmit j m ijt (7.5) 0,0 =?=? tmAVAV m i m it (7.6) 0, =???? tWiLocVPcapAV timitmit (7.7) Equation (7.5) is the conservation of flow for the vehicles and mitAV is the number of vehicles of mode m at node i for time t = 0. Equation (7.6) limits the total number of vehicles of each mode available at time t = 0. Equation (7.7) prevents the temporary facilities that are not open from accepting any vehicular flow and enforces the vehicle parking capacity for the facilities that are open. 7.1.2 Numerical Experiments In the numerical experiments that were presented in previous chapters, it was always assumed that the supplies were stored only at the permanent source nodes (i.e. LC, CSS, VEN). Also the initial distribution of supplies among these 3 sources was a given fixed data. For example for most numerical cases in this study, it was assumed that 40% of supply was stored at the logistics center site in Atlanta, GA. 20% was stored at commercial storage site in Charlotte, NC and 40% of supply was stored at the vendors site in Nashville, TN. A mostly similar approach was used to initially locate the trucks and planes. That prepositioning scheme was arbitrary and obviously can be very far from optimal. 188 In this section a number of cases are generated and solved in order to investigate the properties of an optimal prepositioning scheme. Prepositioning is considered at different levels. First, it is assumed that only the 3 permanent sources can be used for prepositioning and the objective is to find the optimal amounts of supply for each commodity to be stored at these 3 facilities. Secondly, this constraint is removed and it is assumed that all of the facilities in the federal level network can be used for prepositioning. Another important issue is the prepositioning of the vehicles. Numbers and locations of the transportation vehicles at the beginning of the operation can also have a major impact on the efficiency and speed of the operations. If the vehicles are not located optimally at the beginning of the operations, a long period of very valuable time is wasted before the vehicles can arrive at the supply sites, load the relief items and start the delivery process. In fact, the prepositioning of supplies and vehicles are two problems that are related very closely. If supplies are located optimally but vehicles are not readily available then the operations cannot start. On the other hand, prepositioning of vehicles without considering the supply sites is not beneficial either. Consequently, it is very important to consider the prepositioning of supplies and vehicles in conjunction. Optimizing the two problems jointly is the only way to find the best prepositioning scheme. Six cases are generated and solved to optimality to test the effects of prepositioning at different levels. Table 7.1 introduces these cases and presents their optimal objective function value. Case 1 is the base case when no 189 preposition is applied. In the base case, 40% of supplies are located at LC, 20% at CSS and 40% at VEN site. In the base case, 60 trucks are available that are evenly distributed among these 3 source nodes. Also, 20 C130 cargo planes are available that are assumed to be initially located at CSS and VEN facilities, 10 planes at each location. The amount of supply for each commodity and also the number of vehicles of each type is the same for all cases. Facility location problem is considered for all 6 cases with a maximum of 2 MOBs, 2 FOSAs and 4 SSAs. Supplies and vehicles can only be prepositioned at an open temporary facility site and is subject to the capacity constraints of each particular facility. Table 7.1 Introduction of prepositioning case studies CASE Description Objective Value Improvement (%) 1 base case - amounts and locations of supplies and vehicles are fixed at time 0 3.50445 0 2 supplies can be prepositioned optimally at 3 source facilities but vehicles are fixed 3.43035 2.1 3 supplies can be prepositioned optimally at 11 federal facilities but vehicles are fixed 3.3051 5.7 4 amounts and locations of supplies are fixed but vehicles are prepositioned freely 2.51325 28.3 5 Supplies can be prepositioned optimally at 3 source facilities, vehicles are prepositioned freely 2.51325 28.3 6 supplies can be prepositioned optimally at 11 federal facilities, vehicles are prepositioned freely 1.75203 50 The main conclusion from the analysis presented in table 7.1 is that the prepositioning of supplies and vehicles can be very effective. All of the cases with prepositioning show an improvement over the base case scenario. Cases 2 and 3 190 show a modest improvement over the base case. The main reason for that is because the initial distribution of the empty vehicles is not optimal. In these two cases, the vehicles are still located at the source nodes, so there is not much advantage in prepositioning the supplies in other nodes. However, in cases 2 and 3, the supply distribution is adjusted to the initial locations of the vehicles and resulted in a better objective function. In cases 4, 5, and 6 the prepositioning of supplies is combined with optimal distribution of empty vehicles. Since both supplies and vehicles are located in optimal locations, the operations can proceed much faster than the base case scenario. The maximum gain is observed in case 6. In case 6, supplies can be prepositioned in any of the 11 federal level facilities and the vehicles are free available at the same 11 sites. Consequently, the supplies are prepositioned at the nodes closer to the affected states combined with the appropriate number of vehicles available at those sites. This resulted in huge saving of 50% in the objective function compared to the base case. Table 7.2 illustrates the percent of total supply that is prepositioned at each facility type for the 6 cases introduced previously. In case 1, supply assignment is fixed for the source nodes. In case 2, all supply must remain in the source facilities however, it is possible to shift it among the source facilities. As a result, the supplies are moved from LC to CSS and VEN mainly because the majority of empty vehicles are stored at these 2 facilities. In case 3, prepositioning in temporary facilities is allowed. As a result, 30% of supply is moved to 2 FOSA facilities that are closer to the disaster areas. In cases 4 and 5 no supplies are 191 stored at MOB and FOSA facilities because there are no empty vehicles available at those points. Joint prepositioning of vehicles and supplies are allowed in case 6. Therefore a large portion of supplies are prepositioned at MOB and FOSA facilities in case 6. Table 7.2 Percent of total supply prepositioned at each facility type % Supply LC CSS VEN MOBs FOSAs Case 1 40.00 20.00 40.00 - - Case 2 17.92 50.00 32.08 - - Case 3 46.83 23.17 - - 30.00 Case 4 40.00 20.00 40.00 - - Case 5 20.04 50.00 29.96 - - Case 6 1.75 14.17 21.50 32.58 30.00 To analyze the effect of prepositioning on the distribution of empty vehicles refer to Table 7.3. In table 7.3, the number and locations of vehicles for ground and air transportation is shown for the 6 different prepositioning cases. In the first 3 cases the numbers and locations were fixed and only shown for comparison. The general pattern in cases 4, 5 and 6 indicates that the majority of trucks are prepositioned at MOB and FOSA facilities but most of the planes remained at the source nodes. The trucks are moved to MOBs and FOSAs to be closer to the disaster states. Because at the state level, only ground transportation is available and the trucks located at FOSAs can continue to the state level networks and eventually deliver the relief items to the final PODs. On the other hand, at the federal level network, distances among the nodes are much longer 192 which is preferable to be covered by much faster air transportation mode. This pattern in prepositioning of the vehicles is a very important observation that shows the successful collaboration between these 2 transportation modes. Table 7.3 Number of vehicles prepositioned at each facility type No. Of Vehicles Ground Transportation Air Mode LC-CSS-VEN MOBs FOSAs LC-CSS-VEN MOBs CASE 1 60 - - 20 - CASE 2 60 - - 20 - CASE 3 60 - - 20 - CASE 4 3 45 12 17 3 CASE 5 4 39 17 16 4 CASE 6 1 35 24 20 - Overall, it is concluded that prepositioning is very important and it can greatly improve the speed and efficiency of emergency response logistics. It is also shown that prepositioning is especially effective when it is considered in the deployment of relief supplies as well as the transportation capacity both at the same time. If prepositioning of supplies is done without considering the availability of transportation means, the improvements are very marginal in the best case. 7.2 Equity As initially introduced in chapter 3, considering equity and fairness among aid recipients is an important issue. Based on the geographical dispersion of victims and availability of resources over time and space, it is easy to favor the 193 demands of one group of victims over another. Some variations are inevitable; however, the ideal pattern is to distribute the help items evenly and fairly among all victims. The mathematical models and procedures with general objective functions are prone to ignore equity and level of service requirements in order to get a better numerical solution. It is very important to realize the need for procedures and constraints that prevent any sort of discrimination among victims, as much as possible. The equity constraint between populations can be defined over time, and over commodities. It is not appropriate to satisfy all the demands of one group in early stages while the other group of victims does not receive any help until very later times. It is more acceptable to fairly distribute the available relief items among all recipients even though it might not be enough to satisfy all demands. The equity over commodities is also important. For example, it is not acceptable to send all the available water to one group of victims and send all the available meals to another group. In chapter 3, equity constraints were mathematically modeled for the first time. Those equity constraints were tested in the preliminary numerical study in chapter 4. The results indicated that the proposed equity constraints were successful in implementing the required minimum level of service for demand points. However, the equity constraints would increase the size of the model and make it much more difficult to solve. In fact, when all equity constraints were used, the CPLEX solver was not able to find the optimal solution in a reasonable time. In this section, it is tried to analyze the equity constraints proposed in 194 chapter 3 and also introduce a new set of equity constraints with stronger properties. 7.2.1 New Equity Constraints Equations 3.24, 3.25, and 3.26 represented the equity constraints in chapter 3. Equation (3.24) enforces a minimum percentage of total demand for a specific commodity c, to be satisfied by the time period t. Equation (3.25) requires that from all commodities being delivered to node i by time t, at least minb percent to be commodity c. Equation (3.26) ensures that the sum of total commodities delivered to point i, to be more than a minimum percentage of the commodities that are being delivered to all other demand points. tcPODiDem X t c ti t m j cm ttji jim ,,min )( ???? ??? ? ? ? ?? a Repeated (3.24) tcPODiX X c t m j cm ttji t m j cm ttji jim jim ,,min )( )( ??????? ??? ? ?? ? ?? b Repeated (3.25) tPODiX X i c t m j cm ttji c t m j cm ttji jim jim ,min )( )( ???????? ???? ? ?? ? ?? g Repeated (3.26) Equation (3.24) enforces a minimum percentage of total demand for a specific commodity c, to be satisfied by the time period t. It might not be always possible to deliver the required amount by time t; in that case, this constraint makes the optimization problem infeasible. 195 Equation (3.25) requires that from all commodities being delivered to node i by time t, at least minb percent to be commodity c. However, if the demand for a specific commodity is only a small portion of the total demand at that node, then minb cannot be enforced and can cause infeasibility. Equation (3.26) ensures that sum of total commodities delivered at point i to be more than a minimum percentage of all the commodities that are being delivered to all other demand points. Eq. 3.26 might also perform poorly when the amount of demands between different PODs has large fluctuations. In the following, improved equity constraints are proposed to deal with these shortcomings. It is important to notice that the amounts delivered at each node should be proportional to their total demands, in other words: ( ) ( ) ( ) ( ) ( ) ( )nndemand delivered demand delivered demand delivered ??? ... 2 2 1 1 (7.8) This translates to ( ) ( ) ( ) ( ) ( ) ( )n ndemand dunsatisfie demand sunsatisfie demand dunsatisfie ??? ... 2 2 1 1 (7.9) We define a new variable UDRi = relative unsatisfied demand at node i (Total unsatisfied demand at node i divided by total demand at node i). Then we have: nUDRUDRUDR ??? ...21 (7.10) To implement this in mathematical programming language, we should define a tolerance factor ?t and enforce it for each pair of demand nodes. It translates to 196 tNjiUDRUDR ttjti ,, ???? m (7.11) Which is equivalent to tjiUDRUDR ttjtit ,,????? mm and finally we have ?? ??? ???? ???? tNjiUDRUDR tNjiUDRUDR t t i t j t t j t i ,, ,, m m (7.12) Equation 7.12 requires the difference of percentage of unsatisfied demand between each pair of demand nodes to be smaller than a tolerance factor ?t. This new formulation compared to the equity constraints 3.26 is simpler, more compact, and can better describe the concept of equity among demand nodes. Also, the issue of infeasibility caused by previous equations is solved since equation 7.12 considers the relative unsatisfied demand. Equation 7.12 is compact but it can still be improved. If there are N nodes in the network, we need to enforce equation 7.12 for each pair of nodes. In fact we will have N*(N-1)*t equations when the number of nodes is N. In this case, the number of constraints grows very fast with number of demand nodes. To deal with this issue we define 2 new variables (Rmin and Rmax) and reformulate equitation 7.12: Rtmin = Auxiliary variable for minimum relative unsatisfied demand among all nodes at time t Rtmax = Auxiliary variable for maximum relative unsatisfied demand among all nodes at time t ??? ??? ? ?? ?? ?? t tt tt tt RR iRUDR iRUDR i i mminmax max min (7.13) 197 Equations 7.13 find the minimum and maximum relative unsatisfied demands and require their difference to be less than a tolerance parameter ? set by the user. This provides an exclusive control that was not possible before. Also, equation 7.13 is more efficient than equation 7.12. If the number of nodes are N, it only requires (2n+1) constraints to formulate equity compared to N*(N-1) constraints in equation 7.12. This is very important especially when number of nodes is very large. For example in case of 100 nodes, equation 7.12 needs 9900 constraints compared to only 201 constraints if equation 7.13 is being used. Using similar analogy for other equity constraints, new equations are derived. Equations 7.14 to 7.20 present the new equity constraints that would replace equation 3.24, 3.25 and 3.26 in the original formulation. Minimum fill rate: tcViDemUD t c it c it ,,).1( min ???? ?a (7.14) Equity over Commodities: tcViDemrUD t c it t t c it ,,.min ??? ?? (7.15) tcViDemrUD t c it t t c it ,,.max ??? ?? (7.16) trr ttt ??? lminmax (7.17) Equity among different Points of Distribution: tViDemRUD c t c it t c t c it ,.min ??? ???? (7.18) tViDemRUD c t c it t c t c it ,.max ??? ???? (7.19) tRR ttt ??? mminmax (7.20) 198 7.2.2 Numerical Experiment To test the effect of using new equity constraints, 3 new numerical cases are generated and solved. The results are compared to each other as well as to the base case scenario. In these cases, only the equity among points of distribution is considered and it is enforced once at t = 12 hours. Case 1, is the base case scenario without the enforcement of any equity constraints. In Case 2, ? = 0.5 is required that means the difference between highest and lowest relative unsatisfied demand should be less than 50%. In cases 3 and 4 higher equity requirements are enforced by reducing the tolerance to ? = 0.25 and ? = 0.10, respectively. Table 7.4 presents the outcome of the optimal solutions for each case. In Table 7.4, the total demand for each of the 20 POD nodes is shown as well as the total unsatisfied demand at each POD after 12 hours of operation for the base case and other 3 cases. It is very interesting to notice that applying equity constraints in case 2 and case 3 changed the details of the operation however it did not change the objective function value compared to the base case. This indicates that the model was capable of satisfying the new level of service requirements without sacrificing the objective function. For case 4, the objective function is a little higher which is the trade off for satisfying the greater restrictions when ? = 0.1. In table 7.4 the initial demand and unsatisfied demand after 12 hours are given for different cases. From these results, we can calculate the relative unsatisfied demand for each node and each case. The relative unsatisfied demand 199 is equal to unsatisfied demand at each node divided by its initial demand. Table 7.5 reports UDR values. Table 7.4 Effect of equity constraints on demand satisfaction at PODs POD Node Initial Demand Unsatisfied Demand after 12 hours Base case Case 2 Case 3 Case 4 22 13500 6625 13500 8470.588 10125 23 36000 12000 18000 22588.24 27000 24 21000 19250 14333.33 13176.47 15750 25 24000 18000 13500 15058.82 18000 26 37500 13500 23166.67 23529.41 28125 27 42000 24500 28500 26352.94 31500 28 55500 41000 31500 34823.53 41625 29 10500 10000 6000 6588.235 7875 30 37500 33500 31500 32904.41 31875 31 25500 19500 13500 22375 21675 32 12000 11000 12000 10529.41 10200 33 16500 16500 9333.333 12247.79 14025 34 28500 27500 28500 25007.35 24225 35 25500 25500 25500 22375 21675 36 39000 39000 33000 34220.59 33150 37 57000 41125 54500 48808.82 48450 38 10500 10500 10500 9213.235 8925 39 46500 34500 40500 40801.47 39525 40 18000 18000 18000 15794.12 15300 41 43500 43500 38666.67 37722.79 36975 Objective Function (E+7) 3.50445 3.50445 3.50445 3.51866 200 Table 7.5 Relative unsatisfied demand for various equity tolerance levels POD Node Relative Unsatisfied Demand after 12 hours Base case Case 2 Case 3 Case 4 22 0.491 1.000 0.627 0.750 23 0.333 0.500 0.627 0.750 24 0.917 0.683 0.627 0.750 25 0.750 0.563 0.627 0.750 26 0.360 0.618 0.627 0.750 27 0.583 0.679 0.627 0.750 28 0.739 0.568 0.627 0.750 29 0.952 0.571 0.627 0.750 30 0.893 0.840 0.877 0.850 31 0.765 0.529 0.877 0.850 32 0.917 1.000 0.877 0.850 33 1.000 0.566 0.742 0.850 34 0.965 1.000 0.877 0.850 35 1.000 1.000 0.877 0.850 36 1.000 0.846 0.877 0.850 37 0.721 0.956 0.856 0.850 38 1.000 1.000 0.877 0.850 39 0.742 0.871 0.877 0.850 40 1.000 1.000 0.877 0.850 41 1.000 0.889 0.867 0.850 Min 0.333 0.500 0.627 0.750 Max 1.000 1.000 0.877 0.850 ? 0.667 0.50 0.25 0.10 201 Table 7.5 shows that for the base case which had no equity restriction, minimum and maximum relative unsatisfied demands are at 33.3 and 100 percent. This means that for the base case, at least one node has 33% unsatisfied demand while at least one other node has 100% unsatisfied demand after 12 hours of operations. This gap is very large and might not be acceptable from the equity perspective. As it is shown at the bottom of table 7.5, for cases 2, 3, and 4 this discrepancy is lower. In case 2, the differences between the nodes with highest and lowest satisfaction rate are 50%. The same measure is 25% and 10% for cases 3 and 4 respectively. Figure 7.1 better illustrates the outcome of using equity constraints. In figure 7.1 the relative unsatisfied demand is depicted for the four cases described earlier and calculated in table 7.5. It is evident that in the base case, when there are no restrictions, the satisfaction rates for different PODs have large fluctuations. The fluctuations among PODs are reduced when equity constraints are enforced. In fact, for case 4 with ? = 0.1, the fluctuations are very much controlled and the differences are reduced to less than 10%. 202 Figure 7.1 Satisfaction rate variations among points of distribution 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 1.1 2223242526272829303132333435363738394041 Re lati ve UD POD node Number Base Case ?=50% ?=25% ?=10% 203 Chapter 8: Summary and Future Research In today?s society that disasters seem to be striking all corners of the globe, the importance of emergency management is undeniable. Much human loss and unnecessary destruction of infrastructure can be avoided with better planning and foresight. When a disaster strikes, various aid organizations often face significant problems of transporting large amounts of many different commodities including food, clothing, medicine, medical supplies, machinery, and personnel from several points of origin to a number of destinations in the disaster areas. The transportation of supplies and relief personnel must be done quickly and efficiently to maximize the survival rate of the affected population. Federal emergency management agency (FEMA) is the primary organization for preparedness and response to federal level disasters in the United States. FEMA has a very complex supply chain structure to provide the disaster victims with critical items after a disaster which involves multiple organizations and spreads all across the country. Unfortunately, inadequate response to hurricanes Katrina and Rita showed the critical need for better mechanisms in emergency operations. In this research, first FEMA?s supply chain structure is investigated. There are seven main components in FEMA?s supply chain to provide relief commodities for disaster victims that are briefly described here: 1. FEMA Logistics Centers (LC): permanent facilities that receive, store, ship, and recover disaster commodities and equipment. 204 2. Commercial Storage Sites (CSS): permanent facilities that are owned and operated by private industry and store commodities for FEMA. Freezer storage space for ice is an example. 3. Other Federal Agencies Sites (VEN): representing vendors from whom commodities are purchased and managed. Examples are Defense Logistics Agency (DLA) and General Services Administration (GSA). 4. Mobilization (MOB) Centers: temporary federal facilities in theater at which commodities, equipment and personnel can be received and pre- positioned for deployment as required. In MOBs commodities remain under the control of FEMA logistics headquarter and can be deployed to multiple states. MOBs are generally projected to have the capacity to hold 3 days of supply commodities. 5. Federal Operational Staging Areas (FOSAs): temporary facilities at which commodities, equipment and personnel are received and pre-positioned for deployment within one designated state as required. Commodities are usually being supplied from MOB Centers, Logistics Centers or direct shipments from vendors. FOSAs are generally projected to hold 1 to 2 days of commodities. 6. State Staging Areas (SSA): temporary facilities in the affected state at which commodities, equipment and personnel are received and pre- positioned for deployment within that state. Title transfers for delivered federal commodities and cost sharing are initiated in SSAs. 205 7. Point of Distribution (PODs) Sites: temporary local facilities in the disaster area at which commodities are distributed directly to disaster victims. PODs are operated by the affected state. The modeling and optimization techniques established in commercial supply chain management seem to be the most relevant approach to be used in emergency logistics. Some recent studies emphasized that some supply chain concepts share similarities to emergency logistics and therefore tools and methods developed for commercial supply chains can be successfully adapted in emergency response logistics. However, using commercial supply chain techniques in disaster management is still in its infancy. The partial reason is the difference in the strategic goals of commercial supply chain with the goals of disaster response logistics. The main goal in commercial supply chain is to minimize the cost or maximize the profit of the operations. Actions are justified if they increase the profit but are not perused if their cost is more than their profit. However, humanitarian organizations are mostly non-profit organizations with the idea of providing critical services to the public in order to minimize the pain and sufferings of the affected populations. There are not many publications that directly applied network modeling and optimization techniques in disaster response. Among those studies, there is no model that has integrated the interrelated problems of large-scale multicommodity multimodal network flow problem, the vehicle routing problem with split mixed pickup and delivery, and the optimal location finding problem with multiple 206 layers. Also to the best of our knowledge, there is no mathematical model that describes the special structure of FEMA?s supply chain system. The goal of this research is to develop a comprehensive model that describes the integrated logistics operations in response to natural disasters at the operational level. The proposed mathematical model integrates three main components. First, it controls the flow of several relief commodities from sources through the supply chain until they are delivered to the hands of recipients. Second, it considers a large-scale unconventional vehicle routing problem with mixed pickup and delivery schedules for multiple transportation modes. And third, following FEMA?s complex logistics structure, a special facility location problem is considered that involves four layers of temporary facilities at the federal and state levels. Such integrated model provides the opportunity for a centralized operation plan that can effectively eliminate delays and assign the limited resources in a way that is optimal for the entire system. The proposed model considers sending multiple relief commodities (e.g., medicine, water, food, equipment, etc) from a number of sources to several distribution points in the affected areas through a chain structure with some intermediate transfer nodes. The supplies may not be available immediately but arrive over time. It is a difficult task to decide on the right type and quantity of relief items, the sources and destinations of commodities, and also how to dispatch relief items to the recipients in order to minimize the total unsatisfied demand for all disaster victims. 207 Vehicle routing and scheduling during the disaster response is also very important. A large number of vehicles are used in response to large-scale disasters. The proposed model is able to keep track of routings for each individual vehicle. Also, the model provides a detailed schedule for pickup and delivery of relief commodities by each vehicle in each transportation mode. Nonetheless, the vehicle routing in disaster situations are quite different from conventional vehicle routings. The vehicles do not need to form a tour and return to the initial depot, but they can be assigned to a new path at any time. They are expected to perform mixed pickup and delivery of multiple items between different nodes of the network as the supplies and demands arise over time. During the initial response time it is also necessary to set up temporary transfer facilities to receive, arrange, and ship the relief commodities through the distribution network. The proposed model considers optimal selection of several facilities that results in the maximum coverage of the affected areas and the minimum delays for supply delivery operations. Usually the number of these temporary facilities is limited because of the equipment and personnel constraints. Each facility in the model is subject to some capacity constraints. Various capacities are defined for operations such as sending, receiving, and storing commodities. These capacities can be different for each facility and are determined based on the type, size and layout of that facility. Also the availability of personnel and equipment may influence the capacities. The capacity constraints are defined in terms of the weight or volume of the commodities as well as in 208 terms of the number of the vehicles that are sent, received, or parked at each facility during each time period. The other important issue is considering equity and fairness among aid recipients. Based on the geographical dispersion of victims and availability of resources over time and space, it is easy to favor the demands of one group of victims over another. Even though some variations are inevitable, the ideal pattern is to distribute the help items evenly and fairly among the victims. The proposed model recognizes this need and considers a set of constraints that prevent discrimination among victims, as much as possible. A set of preliminary numerical experiments was designed to test the proposed formulation and evaluate the properties of the optimization problem. Different case studies were generated based on the same structure of an imaginary hurricane scenario to analyze the effects of different parameters. In general, the proposed modeling framework produced reasonable outcomes. It was able to provide the level of details required in the disaster response logistics at the operational level. For simple cases and small size problems, the commercial solver was able to find the optimal solutions; however, for the full size problem CPLEX commercial solver was unable to deliver good results within a meaningful computation time. It is concluded that better solution algorithms or heuristics are needed to address the larger problem instances or real world size problems. To develop solution algorithms, first some solution approaches for general integer programming from previous studies in the literature were reviewed. It was 209 concluded that the current model is very complex and a reliable exact solution method is not available that would be computationally attractive or affordable. Consequently, it is more realistic to develop fast and efficient heuristic algorithms to find near optimal solutions. To solve the proposed mathematical model, the problem was structurally decomposed into three sub-problems: multicommodity network flow, location finding with multiple layers, and general integer vehicle routing problem. These 3 problems were solved one-by-one, however the interrelation between these 3 problems were preserved at all times. Multicommodity network flow problem is a linear program and is considered easy since efficient commercial solvers exist that can solve very large LP programs quickly. To solve the multi-layer facility location problem, 4 heuristic methods are proposed. From those, the Branch-and-bound with hierarchical decomposition and the highest capacity ratio were the 2 algorithms that showed better results. To solve the general integer vehicle routing problem, four heuristic algorithms were proposed. The algorithms were tested with large-size numerical experiments. All four algorithms were successful in finding good integer solutions. The convergence rates of the proposed algorithms were also much faster than the commercial solver for the same optimality gap. The proposed VRP algorithms were compared to each other. It was concluded that Y-list and Y-list Modal algorithms were better in solution quality and the convergence speed. However, when worst-case scenario is considered, T- 210 counter and Origin-based T-counter algorithms were shown to have better performances. Finally, all four algorithms were used to solve several random generated problem instances. It was concluded that the proposed algorithms could find good solutions very quickly. In fact, for most cases less than 2% optimality gap was reached in less than 2 minutes of computation time. In chapter 6, in-depth analyses of different aspects of the proposed mathematical model were presented. These analyses were provided in order to better illustrate the capabilities of the model and also examine model?s sensitivity in various circumstances. The analyses were divided into two main categories: sensitivity analysis of the structural parameters of the model, and sensitivity analysis of the main input values of the model. It was shown that changing some input parameters not only affects the optimization results but it can also largely change the problem size and solution computation times. The structural parameters included number of commodities C, number of transportation modes M, and the time-step resolution t. To test the effect of the number of commodities, four cases are considered each with one, two and three and four commodities respectively. It is shown that the number of variables and number of constraints increases linearly with the number of commodities. Number of variables grows faster than number of constraints. More importantly, it is shown that when number of commodities increase, the problem becomes much harder to solve and the computation time increases rapidly. This increase in computation time was exponential in the range tested. 211 Number of transportation modes is also an important factor. Using multimodal compared to single modal transportation not only increases the size of the formulation but also makes the problem much harder to solve. In multimodal systems there is this question of which transportation mode to utilize and how to balance the commodity flows among different modes. Another concern in multimodal environments is considering the intermodal transfer. From application?s perspective it is important to provide suitable facilities and the required equipments in order to transfer the relief commodities between transportation modes quickly and efficiently. From modeling perspective it is important to consider the properties of each transportation mode and correctly model the delays during intermodal transfers. Multimodal problem is more difficult to solve. In fact, for a numerical example the number of variables and constraint in multimodal case is only about 10% more than single-mode case but it was much more difficult to solve the multimodal problem and it took about 120 times longer to solve the multimodal problem. The other structural parameter in the model is the length of time-steps t. Time-step t is the length of time between two consecutive states that the problem is being modeled. Selection of appropriate time-step is a very important factor that can affect the performance and accuracy of time-space networks dramatically. When t is short, the situation on the ground can be modeled in greater details which would not be possible with longer time-steps. So from accuracy perspective, it is favorable to have shorter time-steps. On the other hand, 212 for each time period in the planning horizon, one layer of physical network will be added to the time-space structure of the problem. This makes the problem size grow extremely fast with the number of time-steps in the planning horizon. Finding a reasonable time-step t is an important modeling challenge. In selecting the appropriate time-step one should consider the level of accuracy that is required for that specific application and also the computational power that is available to them. For any specific application, it is recommended to perform an initial analysis and then select the appropriate time-step length based on the required accuracy and availability of computational resources. The three main categories of input parameters are: parameters of the facility location problem, parameters of the vehicle routing problem and parameters of the commodity flow problem. It is shown that the number, location and capacity of the facilities have a major affect on the emergency response operations. In analysis of number of facilities, it was concluded that adding more facilities is beneficial in reducing total unsatisfied demand. However, these improvements became marginal when more that 2 MOBs, 2 FOSAs and 4SSAs are selected. The loading and unloading capacity in each facility is shown to impact the flow of commodities. Numerical studies tested a range of capacities for both loading and unloading capacity factors. It was concluded that investments to expand the capacity should improve both capacities at the same time. If one of the capacities is kept at a fixed level, then additional capacity of the other type remains unused. 213 For vehicle routing problem, number of available vehicles is one of the main factors. It is concluded that the number of vehicles not only affects the model?s results, but it also affects the difficulty of the model and CPU time to find the optimal solution. When the number of vehicles is very low or very high, it is much easier and faster to solve the model. For an in-between range of vehicle numbers, it can become very difficult and time consuming to find the optimal solution. This range is problem-specific and can depend on the other model inputs as well. Researchers and practitioners should be aware of this behavior and perform the similar analysis for a range of vehicle numbers that is appropriate for their specific application. Another factor that can affect the performance of the entire model is the capacity or type of vehicles that are used during the response operations. The general conception is that the higher capacity is always better. However it might not always be possible to use the largest vehicles in the fleet. Analysis of capacity of trucks and planes in the system indicated that the problem was more sensitive to the capacity of trucks than to the capacity of planes. This might be simply due to the fact that the ground transportation does the majority of the. The other factor in vehicle routing problem that affects the response operations is the travel speed. Faster vehicles are favorable from two perspectives. First, the flow of relief commodities through the network can happen faster. Second, empty vehicles can travel faster and reach the pickup nodes to start another round of deliveries in a shorter period of time. However, the travel speed might be reduced in the disaster area due to the inclement weather or road 214 blockings. It was shown that lower travel speeds will increase the unsatisfied demands. Comparing the travel speeds, it was shown that travel speed at federal level network had a stronger effect than travel speed at the state level network. In other words, if we could choose to improve the conditions of the roads in either federal level or state level roadways, it would be more rewarding to improve the federal level links. For commodity flow problem, sensitivity analysis was performed to investigate the effect of various amounts of supply, random demands and the relative urgency factor. In analyzing available supplies, it was shown that the lack of supplies at source nodes had a large negative effect on the objective function value. On the other hand, when supplies were abundant, the objective function was only slightly improved. In fact, other limitations such as transportation capacity and facility constraints limit the amounts of supply that could be delivered and having extra supply is not always beneficial. Locations and amounts of demands are the other factors that affect the details of the response operations. Variability in demand locations and amounts is a negative factor for emergency response operations. Several cases with random demand values were generated and solved to optimality. It was shown that in spite of variation in demands, the proposed model was successful in managing the demand variations. Also, the same set of facility locations were used for all random cases that could be interpreted as a good measure of robustness of the model in case of fluctuations in demand. 215 Finally, the effect of relative urgency factor was tested. It was shown that using relative urgency factor can be very effective when there was a reason to give priority to a commodity over another. The demands for commodities with higher RU factor were satisfied at very earlier times compared to the other commodities. Use of RU factor is highly recommended when there is a small amount of demand for a commodity with high priority (e.g. medical supplies).It is very effective to use a higher RU factor in the model to make sure that the priority is given to that commodity regardless of its small quantity. RU factor can also be applied to give priorities to some demand points over the others. Numerical experiments indicated that RU factors can be effectively used if a some PODs have higher priorities. Recommendations for Future Research In this section some recommendations are listed for future research. 1- Demand Estimation The mathematical model presented in this dissertation is intended for modeling and optimization of disaster response logistics at the operational level. Consequently, the locations, types, and amounts of demand are assumed to be known at the any time. However, in reality the demands are not known in advance. In fact, the demands might be very different in any special case based on the type and intensity of the disaster and the characteristics of the impacted community. It is recommended to research and develop demand estimation models for all potential disaster types and intensities for the targeted communities. Demand 216 estimation models should predict the potential types and locations of demands. It is also very important to predict the variations of demands over time. 2- Modeling intermodal transfer terminals The model proposed in current dissertation is capable of considering multiple transportation modes. It was assumed that relief commodities can be transferred between different modes at some intermediate facilities. Loading, unloading, and storage capacities are considered for these intermodal transfer terminals. Also, the model considers a fixed delay when commodities from one transportation mode are being transferred to another independent of the volume or the type of commodities. Therefore, this is a relatively simplified version of what happens in real world scenarios. It is recommended to do further research on the operational details of intermodal transfer terminals. These terminals are large-scale facilities with various systems and mechanisms inside them. It is suggested to try to model the interactions inside these facilities and then possibly combine the resulted models with the integrated logistic model proposed in this research. 3- Comparison to real world disaster scenarios The integrated mathematical model proposed in this research is a general framework that can be adapted to different disaster scenarios. However, every disaster scenario can be different based on its type, intensity, geographical region and also the amount of available resources and infrastructure. It is recommended to do further research on different disaster types and investigate the requirements specific to that disaster. By doing so, disaster 217 response organizations can improve their operations and gain valuable information in order to deal with future disasters. It is also recommended to gather data and information from previous disasters and use them to calibrate the model proposed in this dissertation in order to tailor this general-purpose model to their specific region and potential disaster scenario. 4- Third party logistics provider and collaborative response In this research, it is assumed that a central organization is in control of managing and operating the entire system. In the United States, federal emergency management agency is the main organization responsible in dealing with large-scale disasters. However in general case, several organizations might run parallel response operations and these organization might act independently. It is also possible that the main relief organization uses other logistics providers (e.g. Contractors, state or local organizations) to help with response to large-scale disasters. It is recommended to investigate the roles and responsibilities of these third party logistic (3PL) providers. It is recommended to develop mathematical models (similar to the one proposed in this dissertation) specific to the operations of 3PL providers. Then, another important research question is how to integrate the operations of the main response organization with those of 3PL providers in order to maximize the benefits for the disaster victims. 5- Dealing with uncertainty 218 In the assumptions section of the proposed model (section 3.4), it was stated that the required data for the model is available and given. For example, it was assumed that following information is given: ? Demands: commodity types, demand locations, demand amounts ? Supply: commodity types, supply locations, supply amounts ? Permanent Facilities: types, locations, capacities ? Temporary Facilities: set of potential sites for each type, capacities of each type ? Network: link-node incidence matrix for each transportation mode ? Vehicles: number of vehicles available for each mode and their initial locations, capacity of each vehicle ? Travel Times: travel time on each link for each transportation mode Since the model is at the operational level, it is assumed that the abovementioned data is given and consequently, the model is deterministic. In the proposed model, the required information is estimated or known at the beginning of the operations and the model can adapt to the new information as the circumstances evolves over time. Therefore, the proposed model in this study is a reactive model and can only adjust to the changes after they happen. An interesting variation of this model is considering a predictive approach. By considering uncertainties and using predicted values for variations of input data over time, the model can plan for events before they happen and achieve greater savings. The investigation and formulation of such a model can be a remarkable contribution. 219 6- Non-homogeneous fleet In the proposed model, it was assumed that all the vehicles of one transportation mode have the same characteristics. The important characteristics are capacity and travel speed of each vehicle type. Modifying the formulation to consider non-homogeneous vehicles within each transportation mode can be investigated in future. Non-homogeneous fleet can be considered in the proposed model however each type of vehicle should be defined as a separate transportation mode. It is suggested to modify the formulation in a way that can consider non- homogeneous fleet of the same transportation mode without defining a new mode that can increase the size of the problem rapidly. 7- Other Objective Functions The objective function modeled in this study was to minimize total unsatisfied demands over time for all commodities and all demand nodes. It is suggested to investigate the possibility of considering and modeling other objective functions. Some examples include maximizing throughput, maximizing utilization, minimizing cost, minimizing operations duration, etc. When such objective functions are formulated, one can analyze and compare the effects of each objective on the details of the operations. 220 Bibliography 1. Alexander, D.E. 1993, ?Natural Disasters?. London: UCL Press 2. Altay, N., and Walter G. G. 2006. ?OR/MS research in disaster operations management?. European Journal of Operational Research, 175, p475-493. 3. Ardekani S. A. and Hobeika A. (1988) ?Logistics problems in the aftermath of the 1985 Mexico City earthquake?. Transportation Quarterly. 42(1), 107-124. 4. Balas, E. ?An Additive Algorithm for Solving Linear Programs with Zero- One Variables?. Operations Research 13(1965):517-546. 5. Barbarosoglu, G. and Arda,Y., ?A two-stage stochastic programming framework for transportation planning in disaster response?. Journal of the Operational Research Society, 2004, 55, 43?53. 6. Barbarosoglu, G., Ozdamar, L. and Cevik, A., ?An interactive approach for hierarchical analysis of helicopter logistics in disaster relief operations?. European Journal of Operational Research, 2002, 140, 118?133. 7. Beale, E.M.L. ?Survey of Integer Programming?. Operation Research Quarterly 16:2(1965):219-228. 8. Beamon, B., 1998. ?Supply Chain Design and Analysis: Models and Methods?, International Journal of Production Economics, Vol. 55, No. 3, pp. 281-294. 9. Beamon, B., ?Measuring supply chain performance?. International Journal of Operations & Production Management; 1999, Vol. 19 Issue 3/4, p275-292. 221 10. Beamon, B. and Kotleba, S., ?Inventory Modeling for complex emergencies in humanitarian relief operations?. International Journal of Logistics. Vol. 9, No. 1, March 2006, 1-18. 11. Beamon, B., ?Humanitarian relief chains: Issues and challenges?, in Proceedings of the 34th International Conference on Computers and Industrial Engineering, San Francisco, CA, 2004. 12. Benders, J.F. ?Partitioning Procedures for Solving Mixed-Variables Programming Problems?. Numerical Methods. 4(1962):239-252. 13. Bodin, L.D., Golden, B.L., Assad, A.A., and Ball, M.O. (1983), ?Routing and scheduling of vehicles and crews. The state of the art?, Computers and Operations Research 10, 69-211. 14. Brown G. G. and Vassiliou A. L. (1993) ?Optimizing disaster relief: real-time operational and tactical decision support?. Naval Research Logistics. 40, l-23. 15. Christofides, N., Mingozzi, A., and Toth, P. (1981), ?Exact algorithms for the vehicle routing problem, based on spanning tree and shortest path relaxations?, Mathematical Programming 20, 255-282. 16. Crainic TG, Rousseau J-M. ?Multicommodity, multimode freight transportation: a general modeling and algorithmic framework for the service network design problem?. Transportation Research B: Methodological 1986; 20:225-242. 17. CRED, ?Annual Disaster Statistical Review: Numbers and Trends 2007?. Center for Research on the Epidemiology of Disasters. http://www.cred.be/ 222 18. Dantzig, G.B., D.R. Fulkerson, and S.M. Johnson. ?Solution of a Large Scale Traveling Salesman Problem?. Operations Research. 2(1954):393-410. 19. Daskin, M.S., 1995. ?Network and Discrete Location: Models, Algorithms, and Applications?. Wiley & Sons, New York. 20. Desrochers M, Lenstra J K, Savelsbergh M W P (1990), ?A classification scheme for vehicle routing and scheduling problems?. European Journal of Operational Research 46: 322?332 21. Dimitruk, Paul. ?Three keys to supply chain management in times of disaster?. Healthcare Purchasing News. Dec 2005. 22. Drezner, Z., Hamacher, H., 2002. ?Facility Location: Applications and Theory?. Springer-Verlag, Berlin. 23. Eksioglu, B., ?Network algorithms for supply chain optimization?. Ph.D. Dissertation, University of Florida 2002. 24. EM-DAT, ?Emergency Disasters Data Base?. http://www.em-dat.net/ 25. FEMA ?Strategic Plan 2008-2013?. http://www.fema.gov/ 26. Fisher, M.L. ?Optimal solution of vehicle routing problems using minimum K-trees?, Operations Research, vol. 42, pp. 626?642, 1994. 27. Fisher, M.L.?The Lagrangian Relaxation Method for Solving Integer Programming Problems?. Management Science 27(1981):1-18. 28. Geoffrion, A.M. ?Lagrangian Relaxation and its Uses in Integer Programming?. Mathematical Programming Study, 2(1974):82-114. 29. Golden, B., L. Bodin, T. Doyle and W. Stewart. ?Approximate Traveling Salesman Algorithms?. Operations Research. 28(1980):694-711. 223 30. Golden, B.L., and Assad, A.A. (1988), ?Vehicle Routing: Methods and Studies?, North-Holland, Amsterdam. 31. Gomory, R.E.?An Algorithm for Integer Solutions to Linear Programs?. In Recent Advances in Mathematical Programming. (R.L. Graves and P. Wolfe, eds.). McGraw-Hill, New York, 1963:269,302. 32. Guelat J., Florian M. and Crainic T. (1990), ?A multimode multiproduct network assignment model for strategic planning of freight flows?. Transportation Science. 24, 25-39. 33. Haddow. G., Jane A. Bullock. J, Coppola. D., ?Introduction to Emergency Management?, Published by Butterworth-Heinemann, 2007 34. Haghani, A. and Oh, S.C., ?Formulation and solution of a multi-commodity, multi-modal network flow model for disaster relief operations?. Transport. Res. Part A-Policy and Practice, 1996, 30(3), 231?250. 35. Harland, C.M., 1996. ?Supply chain management: relationships, chains and networks?. British Academy of Management 7 (Special Issue), S63-S80. 36. Hughes, M.A., 1991. ?A selected annotated bibliography of social science research on planning for and responding to hazardous material disasters?. Journal of Hazardous Materials 27, 91?109. 37. Kemball-Cook, D., Stephenson, R. (1984), ?Lessons in logistics from Somalia?, Disasters, Vol. 8 pp.57-66. 38. Knott R. (1987), ?The logistics of bulk relief supplies?. Disasters 11, 113-115. 39. Knott R. (1988), ?Vehicle scheduling for emergency relief management: a knowledge-based approach?. Disasters 12, 285-293. 224 40. Land, A.H. and A.G. Doig. ?An Automatic Method for Solving Discrete Programming Problems?. Econometrica 28(1960):497-520. 41. Laporte, G., and Nobert, Y. (1987), ?Exact algorithms for the vehicle routing problem?, in: S. Martello, G. Laporte, M. Minoux and C. Ribeiro (eds.), Surveys in Combinatorial Optimization, North-Holland, Amsterdam, 147-184. 42. Laporte, G., ?The vehicle routing problem: an overview of exact and approximate algorithms?, European Journal of Operational Research, vol. 59, pp. 345?358, 1992. 43. G. Laporte, ?Vehicle routing?. In Dell?Amico, Maffioli & Martello (Eds.) Annotated Bibliographies in Combinatorial Optimization. New York, Wiley, 1997. 44. Magnanti, TL, 1981. ?Combinatorial Optimization and Vehicle Fleet Planning: Perspectives and Prospects?. Networks 11, 179-2 14. 45. Mentzer, J.T., W. DeWitt, J.S. Keebler, S. Min, N.W. Nix, C.D. Smith, Z.G. Zacharia. ?Defining Supply Chain Management?, Journal of Business Logistics, (22:2), 2001, pp. 1?25. 46. Nemhauser G.L. and Wolsey L.A. (1999). ?Integer and Combinatorial Optimization?, John Wiley, New York. 47. New, S.J., 1997. ?The scope of supply chain management research?. Supply Chain Management 2 (1), 15-22. 48. New, S.J., Payne, P., 1995. ?Research frameworks in logistics: three models, seven dinners and a survey?. International Journal of Physical Distribution and Logistics Management 25 (10), 60-77. 225 49. Oh, S.C. and Haghani, A., ?Testing and evaluation of a multi-commodity multi-modal network flow model for disaster relief management?. J. Adv. Transport., 1997, 31(3), 249?282. 50. Oloruntoba, R. and Gray, R., ?Humanitarian aid: an agile supply chain?, Supply Chain Management, 2006, 11(2), 115?120. 51. Owen, S. H. and M. S. Daskin, 1998, ?Strategic Facility Location: A Review?, European Journal of Operational Research, 111, 423-447 52. Ozdamar, L., Ekinci, E. and Kucukyazici, B., ?Emergency logistics planning in natural disasters?. Annals of Operations Research, 2004, 129, 217?245. 53. Pan American Health Organization. ?Natural disasters-protecting the public?s health?. Washington DC. (2001) 54. Papadimitrou, C.H. and K. Steiglitz. ?Combinatorial Optimization: Algorithms and Complexity?. Prentice-Hall, Englewood Cliffs, NJ, 1982. 55. Ray J. (1987), ?A multi-period linear programming model for optimally scheduling the distribution of food-aid in West Africa?. MS-Thesis, University of Tennessee, Knoxville, TN. 56. Rawls, C.G., Turnquist, M.A., 2010, ?Pre-positioning of emergency supplies for disaster response?, Transportation Research Part B, Vol44, Issue 4, 521- 534 57. ReVelle, C. S., Eiselt, H.A., and Daskin, M. S., 2008, ?A Bibliography for Some Fundamental Problem Categories in Discrete Location Science,? European Journal of Operational Research, 184:3, 817-848. 226 58. Schilling, D.A., Jayaraman, V. and Barkhi, R., 1993, ?A review of covering problems in facility location?, Location Science 1 (1) (1993) pp. 25-55. 59. Scott, C., Westbrook, R., 1991. ?New strategic tools for supply chain management?. International Journal of Physical Distribution and Logistics 21 (1), 23-33. 60. Tan, K.C., 2001. ?A framework of supply chain management literature?. European Journal of Purchasing & Supply Management 7, 39?48. 61. Thomas, A.S., ?Humanitarian logistics: enabling disaster response?, Fritz Institute, 2007. 62. Thomas, A.S. and Kopczak. L.R., ?From logistics to supply chain management: The path forward in the humanitarian sector?. Fritz Institute, 2005. 63. P. Toth, D. Vigo (Eds.) (2002). ?The Vehicle Routing Problem?, SIAM Monographs on Discrete Mathematics and Applications, Philadelphia. 64. Van Wassenhove, L.N., ?Humanitarian aid logistics: supply chain management in high gear?, Journal of Operations Research Society, 2006,57(5), 475?489. 65. Yi. W. and Ozdamar, L., ?A dynamic logistics coordination model for evacuation and support in disaster response activities?. European Journal of Operations Research 179 (2007) 1177-1193. 66. Zanakis, S.H. and J.R. Evans. ?Heuristic Optimization: Why, When and How to Use it?. Interfaces. 11(1981):84-90.