# LIST OF TABLES viii

THE EXPLOSIVE SCANNING DEVICES ALLOCATION PROBLEM FOR AIRPORT SECURITY SYSTEMS by JACKRAPONG ATTAGARA, M.S. A DISSERTATION IN INDUSTRIAL ENGINEERING Submitted to the Graduate Faculty of Texas Tech University in Partial Fulfillment of the Requirements for the Degree of DOCTOR OF PHILOSOPHY Approved John E. Kobza Chairperson of the Committee Milton L. Smith Co-Chairperson of the Committee Surya D. Liman Christopher J. Monico Sheldon H. Jacobson

Accepted John Borrelli Dean of the Graduate School May, 2006

Copyright 2006, Jackrapong Attagara

ACKNOWLEDGMENTS

I first express my gratitude to my advisor, Dr. John E. Kobza, for his guidance, encouragement, and financial support throughout the whole process of conducting this research. He has been very helpful since the initiation of this dissertation, and has patiently listened to and promptly assisted my problems throughout the dissertation. I am very thankful to have been able to work with him. I thank Dr. Milton Smith for his countless suggestions during my study at Texas Tech University including being cochairperson of my dissertation. I thank Dr. Surya D. Liman, for his helpful advice and for serving as one of my committee members. I am also indebted to, Dr. Christopher J. Monico, for his thoughtfully assistance, especially his suggestions for the construction of Theorem 3.1. I particularly appreciate Dr. Sheldon H. Jacobson, who has served on my committee, and provided many useful and essential suggestions in developing this research. I also wish to thank Dr. Simon M. Hsiang for his valuable advice, and serving as my committee during on qualifying exam. Many people have supported me during my study. In particular, I am deeply indebted to my beloved parents, Karnchit and Viboonporn, my two brothers, Pornpoj and Peera, my only sister, Piyaphan Attagara, for their unfailing love, support, and encouragement during my study abroad. Lastly, special thanks for the inspiration and financial support of my aunt, Dr. Kingkeo Attagara. Without her support, I would not have been able to complete my doctoral studies. I dedicate this dissertation to them.

ii

Last but not least, I acknowledge with gratitude the financial support of National Science Foundation as this research has been supported by the National Science Foundation (DMI-0114046).

iii

TABLE OF CONTENTS

ACKNOWLEDGMENTS ABSTRACT LIST OF TABLES LIST OF FIGURES CHAPTER I. INTRODUCTION 1.1 Motivation 1.2 History of Aviation 1.3 Passenger Screening 1.4 Problem Statement 1.5 Dissertation Outline II. LITERATURE SURVEY 2.1 Introduction 2.2 Aviation Security 2.3 Complexity Theory 2.4 Optimization Problems and Algorithms 2.4.1 Integer Nonlinear Programming Problems 2.4.2 Knapsack Problems 2.4.3 Reliability Problems III. EXPLOSIVE SCREENING DEVICE ALLOCATION MODEL FOR A SINGLE AIRPORT 3.1 Introduction 3.2 Single Airport Scenario 3.3 Notation 3.4 Single Airport Model 3.5 Computational Complexity of the Model

ii vii viii x

1 1 1 4 7 10 11 11 11 17 19 20 22 23

26 26 27 31 32 38

iv

IV.

METHODOLOGY AND EXAMPLES FOR SINGLE AIRPORT MODEL 4.1 Introduction 4.2 The Partial Enumeration Method 4.2.1 Partial Ordering and Numerical Ordering of Binary Vectors 4.2.2 Skipping Rules 4.3 Model Preparation 4.4 Data Generation 4.5 Computational Examples

47 47 48 49 50 53 58 60 66 66 67 68 69 70 74 74 75 76 78 84 90

V.

MULTIPLE AIRPORT PROBLEM 5.1 Introduction 5.2 Multiple Airport Scenario 5.3 Multiple Airport Model 5.3.1 Notation 5.3.2 The Model

VI.

MULTIPLE AIRPORT METHODOLOGY 6.1 Introduction 6.2 The Effect of Interaction Variables 6.3 Airport Greedy Approach 6.4 Variable Greedy Approach 6.5 3-Phase Approach 6.6 Chapter Summary

VII.

EXPERIMENTAL RESULTS FOR MULTPLE AIRPORT PROBLEM 7.1 Introduction 7.2 Exact Algorithm 7.3 Upper Bound 7.4 Experimental Result 7.4.1 Small Size Problem 7.4.2 Large Size Problem 7.5 Chapter Summary

91 91 91 93 97 98 100 129 130 130 132

VIII.

CONCLUSIONS AND FUTURE RESEARCH 8.1 Conclusions 8.2 Directions of Future Research

v

REFERENCES APPENDIX A. MATLAB SOURCE CODE B. EXAPLE OF AIRPORT DATA C. SAS SOURCE CODE

134

138 163 167

vi

ABSTRACT

This dissertation investigates passenger screening at the airport. Before boarding an aircraft, passengers are divided into many classes by a passenger prescreening system according to their perceived risk levels. Subsequently, a passenger is screened by one or more procedures. This dissertation introduces an explosive scanning device allocation model for each group of passengers. The optimization model for multiple-airport security systems is developed, where the focus is on carry-on baggage screening. The objective of the model is to assign both the type and number of devices to each group so that the total security is maximized subject to budget, resource and throughput constraints. This research proves that the model is NP-hard, and introduces heuristics to obtain feasible solutions. Numerical results of the example are provided.

vii

LIST OF TABLES

3.1 4.1 4.2 4.3 4.4 7.1 7.2 7.3 7.4 7.5 7.6

Model Comparison Screening Devices Data Passenger Classes Data Computational Results Comparison of Two Algorithms Computational Time of Lawler and Belle’s Algorithm Three Screening Device Types Data Airports Data Comparison of Upper Bounds and Optimal Solutions for 5 Airports Percentage Error of Upper Bound Computational Time and Objective Value of Four Methods for Five-Airport Problem Percentage Error of Four Heuristics for Five-Airport Problem Data for Seven Types of Screening Device Comparison of Four Heuristics versus Number of Various New Devices for 100-Airport Problem. Percentage Deviation from Upper Bound versus Percentage of the Total Demand Total Demand for Three Test Problems Computational Time of Four Heuristic versus the Number of Device Types for 100-Airport Problem Objective Value of Four Heuristics and the Upper Bound versus the Number of Device Types for 100-Airport Problem

30 60 61 63 64 93 94 94 95 96

98 98 100

7.7 7.8 7.9

101

7.10

103 107

7.11 7.12

108

7.13

109

viii

7.14

Percentage Deviation from Upper Bound of Four Heuristics versus the Number of Device Types for 100-Airport Problem Percentage Deviation for Thirty Data Set of Four Heuristic with Number of Available Devices = 65% of the Total Demand Statistic Data of the Experimental Result in Table 7.15 Shapiro-Wilk P Value for Normality Assumption of Four Heuristics The 95% Confidence of the Four Heuristics’ Percentage Deviation Experimental Result of Four Heuristics when the Percentage of Total Demand = 85 Best Heuristic in Various Conditions Example of Randomly Generated Data for 100 Airports

109

7.15

111 112 113 113

7.16 7.17 7.18 7.19

118 128 164

7.20 B.1

ix

LIST OF FIGURES

3.1 4.1 6.1 6.2 6.3 7.1 7.2 7.3

Example of a Combination of Three Device Types Flowchart of Lawler and Belle’s Algorithm Algorithm AGH Control Abstraction Algorithm VGH Control Abstraction Algorithm 3PH Control Abstraction Percentage Error of Upper Bound Percentage Error of Four Heuristics for the Five-Airport Problem Percentage Deviation from Upper Bound versus Number of New Device Available Upper Bound versus Number of New Device Available Percentage Deviation from Upper Bound of Four Heuristics versus the Number of Device Types Histogram of Residuals for Single Factor Model Residuals versus Heuristic for Single Factor Model ANOVA Result for Comparing VGH, 3PH-A, and 3PH-V SAS Result for Paired Data Histogram of Residuals for Three-Factor Model Residuals versus Factor ANOVA Result for Three-Factor Model from GLM Procedure

35 62 77 83 89 97 99

104 105

7.4 7.5

110 115 117 118 120 123 124 127

7.6 7.7 7.8 7.9 7.10 7.11 7.12

x

CHAPTER I INTRODUCTION

1.1 Motivation Every airplane traveler wants safe and convenient transportation. An aircraft is vulnerable to both unintentional danger, such as a crash due to an engine failure, and intentional danger, such as bombing or a terrorist attack. To prevent the latter, many attempts have been made to improve the level of aviation security. Examples include developing a screening technology to improve the ability to identify potential threats, increasing the number of screening devices and screeners to handle more passengers, and issuing new regulations such as positive passenger bag match. As a result, aviation security relies heavily on screening passengers and baggage before they enter the aircraft. This research examines this extremely important area of aviation security.

1.2 History of Aviation In November 1909, the first passenger airline, Deutsche Luftschiffahrt Aktien Gesellschaft (DELAG), was founded by a German Cavalry General, Count Ferdinand von Zeppelin and a group of entrepreneurs. DELAG used airships to transport roughly 35,000 passengers from 1909 to 1914 (Royal Air Force Museum, 2002). The first regularly scheduled passenger airline using airplanes was established by the Beonist Company, in March 1914 (Cobb and Primo, 2003). The flight route was between Saint Petersburg and Tampa, Florida under the name St. Petersburg-Tampa Air Line (also

1

called the "Airboat Line"). Even though the airline operated only three or four months, the two roundtrips flights per day significantly inspired commercial aviation. It was not until 1931 that the first aircraft hijacking was recorded; a group of armed revolutionaries took control of a Pan American mail plane (U.S. Centennial of Flight Commission, 2004). The first recorded bombing was in November 1955, when Jack G. Graham placed a bomb in his mother’s baggage resulting in the death of all 44 passengers (U.S. Centennial of Flight Commission, 2004). A similar tragedy occurred in January 1960, when a National Airlines DC-6B crashed as a result of dynamite explosion (U.S. Centennial of Flight Commission, 2004). In both incidents, the intent was to kill a specific passenger and collect on their life insurance policy. These bombings sparked demands to limit the maximum amount of life insurance collectable for death in an aircraft accident and to use devices to screen bags for explosives. The number of hijackings increased dramatically during the late 1960s and early 1970s. As a result, President Richard Nixon established an Anti-Hijacking Program to address hijacking in the U.S. The program initiated the use of Federal Air Marshals (FAMS), specially trained marshals used to ensure security during the flight. It also required air carriers to screen all passengers using at least one of the following methods: hijacker profile, metal detector, identification check and physical search. Consequently, the number of hijackings decreased (National Research Council, 1996). In December 1988, a bomb destroyed Pan Am Flight 103 over Lockerbie, Scotland, killing 259 people on board and 11 on the ground. Investigation showed that a bomb made from a plastic substance, was loaded onto the plane in Frankfurt, Germany.

2

The incident shifted airport security attention from firearms to explosive threats. In 1990, the Federal Aviation Administration (FAA) was required to research and develop an explosives detection system (EDS) in order to prevent explosives being loaded onto an aircraft (U.S. Aviation Subcommittee, 2001). During the 10 years following 1991, there was not a successful hijacking in the U.S. This peaceful period was destroyed by a group of terrorists on September 11, 2001. These terrorists hijacked four aircraft to use as weapons to attack the World Trade Center twin towers and the Pentagon, causing the death of thousands. In November 2001, President George W. Bush signed the Aviation and Transportation Security Act (ATSA) in response to the terrorist attacks. The Act included the creation of the Transportation Security Administration (TSA) under the Department of Transportation (DOT). In March 2003, TSA was transferred to the Department of Homeland Security. TSA’s main mission is “to protect our nation’s transportation systems –aviation, waterways, rails, highways, public transit and pipelines – to ensure freedom of movement for people and commerce.” Before September 11, only a subset of checked baggage was screened by EDS before being loaded onto an aircraft, and by 2013, the FAA planned to deploy the federally certified EDS machines to ensure 100 percent screening of checked baggage. Under ATSA, the TSA became responsible for deploying screening machines, and the deadline was shortened from 2013 to 2002. Moreover, TSA was tasked with many major security duties, including developing a basic screener training program, upgrading existing screening technologies, deploying other innovative systems, and overseeing the

3

federalized screener workforce. These security changes required short-term and longterm action. TSA is also developing new technologies and procedures that may be implemented in the future.

1.3 Passenger Screening The main purpose of passenger screening is to prevent a threat from reaching an aircraft. A threat is defined as a danger to an aircraft, including explosive objects, weapons, or terrorists. There are three ways that a passenger can bring a threat onto an aircraft: in their checked baggage, in carry-on baggage, and on their person. In order to detect a threat, all three possibilities must be checked. For checked baggage, the primary procedures used are: explosive detection systems (EDSs), explosive trace detection (ETD), and hand search. An EDS is a very large, sophisticated, and expensive screening machine that uses computerized tomography (CT) to generate three-dimensional images of a bag’s contents. However, EDS is not a flawless machine, its false-alarm rate is somewhat high and the throughput is low (Butler and Poole, 2002). ETD is used in airports where the volume of passengers does not justify the use of EDSs. The ETD operator wipes the bag to pick up samples of chemical from its surface, then inserts the sample into the ETD machine, which looks for the presence of dangerous materials. There are three ETD methods used to sample checked baggage: closed-bag trace, open-bag trace, and non-directed trace. The first method samples only from the outside of the bag and takes an average 47 seconds per bag. The second method searches the inside lining and the outside of the bag and takes 2

4

to 2.5 minutes. The third method wipes every item in the bag larger than a soft-drink can and requires 3 to 4 minutes. The drawbacks of ETD are that it requires more screeners than EDS and the processing rate is quite low. Another security procedure is Positive Passenger Bag Match (PPBM). All U.S. carriers are required to verify that checked baggage is transported with the owner on international flights, including the connecting flights. The checked baggage will be removed from the aircraft if the passenger has not boarded on that flight. This can prevent a terrorist from placing a bomb onto the aircraft in their checked baggage and then leaving the aircraft, but it offers no protection against the suicidal terrorist. Automated Xray and backscatter X-ray machines are additional method for screening checked baggage, but are not certified by TSA. The terrorist can also carry a threat in their carry-on baggage. Generally, carry-on bags are placed on a conveyer belt and inspected by an x-ray machine. An operator then looks for suspicious items. Additional search such as detailed hand search and open bag trace are performed when necessary. Passenger inspection can be classified into two types: passenger profiling or passenger prescreening and physical screening. The FAA developed and implemented the first Computer Assisted Passenger Prescreening System (CAPPS) in 1997. It evaluates passenger information and determines whether a passenger is known or suspected to pose a threat to an aircraft. After September 11, 2001, TSA was required to develop a successor system of CAPPS, designated as CAPPS II. This new version required both passenger details (such as phone number, car ownership, and magazine subscriptions)

5

and minute details (seating records of almost every US passenger), and then assigned a threat assessment value for each passenger according to the level of threat posed to an airplane (Barnett, 2004). CAPPS II was criticized by Privacy Advocates and Congress for being too invasive in the passenger data needed. The TSA officially ended CAPPS II and replaced it with a less intrusive system called “Secure Flight”. The new system will require only passenger name record (PNR) data from air carriers when a passenger makes a reservation. The PNR data varies among air carriers, but usually includes (1) passenger name; (2) reservation date; (3) travel agency; (4) travel itinerary information; (5) method of payment; (6) flight number; and (7) seating number. The system will use an evaluation algorithm to run the PNR against the watch list of the TSA. The list includes individuals labeled as “terrorist” or “no-fly” and “automatic selectee” from the Terrorist Screening Database (TSDB) maintained by the Terrorist Screening Center (TSC). A selectee can be defined as a passenger who is likely to pose a threat and requires special attention. Secure Flight also includes two important features: a redress mechanism and a random element. The redress mechanism will allow selectee passengers to resolve questions if they believe they are unfairly selected for extra security. For the random element, Secure Flight will randomly select additional passengers for additional screening in order to protect against those who try to use the system algorithm to avoid screening. When fully operational, the Secure Flight system is expected to cut down the selectee rate from 15% to 5% or 6%, consequently speeding up the screening procedures at the airport (Singel, 2004).

6

Another mode of passenger inspection is physical screening. Passengers are required to walk through a portal to detect metallic objects. Further search also can be conducted, especially, if there is an alarm from the portal. The process includes hand wand metal detectors, a physical pat down, and x-ray imaging.

1.4 Problem Statement Many changes have been made to improve aviation security and to restore public confidence. These include major policy changes and new inspection technology. Nevertheless, there are many limitations to improving aviation security. Implementing an aviation security system requires a heavy monetary burden. For example, purchasing a single EDS machine costs nearly $2 million, not including facility modification. Over 6,000 EDSs were deployed at a cost of over $12 billion to meet the 100% checked baggage screening requirement. Even if budget were not an issue, it is still impossible to use EDSs in all airports because of its size and weight. Passenger prescreening is inexpensive compared to a physical inspection, but there is no hope that it can reach 100% percent accuracy in the near future. In fact, no single security screening system offers the desired level of efficiency and certainty. A combination of physical inspection and passenger prescreening is required to reach a preferred security level. Aviation security is an ever-changing system. Terrorists are always looking for new methods to attack the current security system. Once the new threats are presented, researchers have to find ways to close the holes in the system. Many technologies are currently under development and likely to be available in the near future to handle newly

7

emerging threats. Therefore, security systems at all the airports in the nation will need renovation from time to time. Although increased security at airports is desirable, it can lead to passenger inconvenience. If passengers experience enough inconvenience due to excessive delay for screening procedures, or intrusiveness, they might change to other modes of transportation or not travel at all. All passengers are required to pass through a series of inspection processes. The time required for each process and waiting time can be significant, especially when the queue grows long. Waiting time in the queue can be reduced by employing an appropriate amount of security equipment and personnel. Security systems should be carefully chosen so that the delay of both selectee and nonselectee passengers is low, and their privacy is not being invaded, while maintaining appropriate security. Over 700 million passengers travel by aircraft in the United States, and more than 1 billion bags are screened each year (Mead, 2002). Meeting the TSA 100% checked baggage screening requirement for 440 commercial airports is very challenging (U.S. General Accounting Office, 2003). Researchers have to determine how many and which type of devices are needed for a particular airport. Furthermore, when a new system is to be deployed, one of the most critical decisions is to determine which airports will receive the system first, and what to do with the existing systems. Several efforts have been made to analyze the performance of aviation security systems. Performance can be measured in many ways, for example, the probability of

8

threat detection (McLay et al., 2006), throughput of the system (Jacobson et al., 2003), and cost effectiveness (Virta et al., 2003). Suppose passenger prescreening such as Secure Flight has been implemented, passengers will be classified into classes. Each class of passengers will receive a different degree of inspection. It should be noted that there are three ways that passengers can carry threats onto aircraft, on their person, in their carry-on baggage or their checked baggage. The explosive screening device allocation problem also can be separated into three independent sub-problems; one for each of the three ways. In this research, we focus on the problem of allocating devices to screen carry-on. The primary objective is to study the explosive screening device allocation problem of multiple airports and present an optimization model to maximize the total security level under a set of realistic constraints. The objective function is derived from the probability of threat detection, while the constraints account for the total budget, system throughput, and resource availability. We first optimize a single airport problem for carry-on baggage, and then expand to cover a group of airports. Initially, the model finds what combinations of the screening devices for each class of passenger work best for, and are appropriate to, an individual airport without considering resource availability. Next, the model is expanded to include a set of airports under consideration instead of just one airport, in one optimization problem. Resource availability is added to the model to determine which airport should receive additional resources. The model is expected to fully utilize existing devices and add new devices only if they are needed to maximize the total security.

9

1.5 Dissertation Outline The remainder of this dissertation is organized as follows. In Chapter 2, the literature related to the proposed research is reviewed and presented in three sections: Aviation Security, Complexity Theory, and Optimization Problems that resemble the Allocation Problem. The development of the model for a single airport is described in detail in Chapter 3. Chapter 4 shows the methodology and results of the experiment for single-airport model. Model for multiple airports is developed in Chapter 5. In Chapter 6, the methodologies for multiple-airport model are presented. Chapter 7 shows the results of the experiment for multiple-airport model. Lastly, the conclusions and future research directions are provided in Chapter 8.

10

CHAPTER II LITERATURE SURVEY

2.1 Introduction Aviation security systems can be improved by several approaches, e.g., upgrading screening device technology and issuing new policies. Many researchers have focused on the deployment of screening devices and the classification of passengers. This research is concerned with the assignment of numbers and types of explosive screening devices to an airport to maximize the total security. The organization of this chapter is as follows. First, an overview of problems in aviation security system is presented. Second, we review complexity theory. Finally, similar types of optimization problems are discussed.

2.2 Aviation Security Following the events of September 11, 2001, the US Government invested heavily in research on screening device technology. No single explosives detection technology can provide the preferred result with confidence and efficiency. Butler and Poole (2002) support the idea of using a passenger prescreening system to distinguish high-risk passengers from other passengers and use EDS (explosives detection system) or EDT (explosive trace detection) for screening checked-baggage of this high-risk group. They point out that using EDS systems alone to screen 100% of the checked baggage is not cost effective. EDS has many flaws, it is large in size, very expensive, and its error rate is nearly 30 percent. The ETD system is very slow and more labor intensive.

11

Moreover, they mention that the current generation of EDS machines is likely to be replaced by a new and better technology for baggage screening in the near future. Poole and Passantino (2003) go into details of the prescreening system suggested by Butler and Poole (2002). They observe that current airport security systems spend too much attention on low-risk passengers and too little on high-risk passengers. They suggest airport security systems should follow a risk-based approach, which means shifting attention from finding threat objects to identifying who is a threat. One alternative of the risk-based approach is to divide all of the airline passengers into three groups as follows, 1. Those about whom a great deal is known, all of it consistent with that person not being a threat. This group includes low-risk passengers, trusted travelers, airline flight crews, especially cockpit crews. (Trusted travelers or registered travelers are people who become members of the a Registered Traveler Program by disclosing their various personal and possibly financial details to the government, These people are pre-cleared and can bypass some inspection stations.) 2. Those about whom very little is known, and what is known is correlated with known risk factors. (Medium-risk passengers) 3. Everyone else. (High-risk passengers) Each group will go through a different level of scrutiny. Passengers in the high-risk group should be inspected by the most sensitive and most costly (probably slowest) explosive detection technologies. Medium-risk passengers are screened at a lower security level than high-risk passengers and low-risk passengers receive the least and fastest attention.

12

In 1997, the FAA developed and implemented the first passenger prescreening system called CAPPS, Computer Assisted Passenger Prescreening System. Its second generation, called CAPPS II, began development after the incident of September 11, 2001. Due to the privacy concern, CAPPS II was replaced before being implemented by a less intrusive system named “Secure Flight” in 2004 (Singel, 2004). Secure Flight analyzes Passenger Name Record (PNR) data and classifies passengers into two groups (U.S. Department of Homeland Security, 2004). The PNR data usually includes (1) passenger name; (2) reservation date; (3) travel agency; (4) travel itinerary information; (5) method of payment; (6) flight number; and (7) seating number. Those who are suspected or known terrorists, and are likely to pose a threat to a system will be put in the ‘selectee’ group. Passengers who are cleared by Secure Flight will be in the second group, namely the ‘non-selectee’ group. The effectiveness of CAPPS II is investigated by Barnett (2004). CAPPS II analyzes the information about a passenger including minute details and calculates a threat assessment for each passenger. He concludes that CAPPS II can work efficiently under certain circumstances but it should not be use as the centerpiece of a security strategy. Determining the number of selectee passengers is necessary for assigning baggage screening security devices in each airport station. Virta et al. (2002) suggest using a mathematical model for calculating the selectee rates for a collection of hub and non-hub airports. They define the selectee rate as the number of selectee passengers divided by the total number of passengers on that flight. The rate depends on the flight

13

origination and destination cities, and characteristics of the passengers, i.e. originating or transferring passengers, and the transfer rate on that flight. Kobza and Jacobson (1997) analyze the performance of security system architectures for both single-device and multiple-device systems. They develop a probability model based on Type I (a false alarm is given) and Type II (a false clear is given) errors. For the multiple-device systems, they consider parallel, series and cascaded system architectures. Their study shows that multiple-device systems outperform singledevice systems for certain error probability measures. This paper is useful in evaluating and comparing the performance of systems, or combinations of devices. When the screening device is used only for selectee passengers, usually, there will be unused capacity. In response to the government mandated policy of 100% screening of checked baggage, Virta et al (2003) use a cost model to study the trade-offs between screening only selectee checked baggage and screening selectee and all or some portion of non-selectee baggage. The associated costs include the cost of detecting and not detecting a threat, fixed, operating and maintenance costs of the machine. Their analysis shows that as the amount of excess capacity used to screen non-selectee baggage increases, the expected annual cost increases, the expected annual cost per checked bag decreases, and the expected annual cost per expected threat detected increases. Continuing the work of Virta et al. (2003) and Kobza and Jacobson (1997), Jacobson et al. (2003) study the cost effectiveness of different configurations of explosive detection devices. They calculate expected annual cost per bag and expected successful threats for different scenarios. Ritchie et al. (2004) extend the work of Jacobson et al.

14

(2003) by applying their model to analyze the cost and expected number of successful threats of various scenarios and the scenario from Poole and Passantino (2003). Rausch (2002) proposes a model where the decision variables are the number of screening devices to install by airport station and time period. Major constraints of the model include demand requirements, and the production schedule for devices. Initially, the objective function maximized the number of bags screened at each station. However, the results exhibited flailing, i.e. a dramatic change of processing rate from month to month, which is not desirable in a real system. An alternative objective function, which minimizes the change from previous plans, is used instead. Lazar Babu (2003) develops a method for grouping passengers and determining the numbers of devices at each check station while not violating the false clears criteria set by TSA. In their problem, if there are m special devices to choose from, for a given n groups of passengers, where 1 ≤ n ≤ m, then there are

? m? ?2 ? ? ? ? n ?

ways of forming groups. Each

group of passengers will have a unique assignment of screening devices. The model decides which way of forming groups to use and determines the percentage of all passengers for each group. The passengers are assigned to each group randomly according to the percentage value. In their model, passengers are not distinguishable. Their objective is to minimize the probability of a false alarm in the airport security system and to prevent screening of all passengers by all check stations. To maximize the total security, McLay et al. (2006) introduce the Multilevel Allocation Problem (MAP), which assigns airline passengers into one of several groups. Each of which has a different configuration of security screening devices. The checked 15

baggage of each group of passengers is then screened by one of the different predefined configurations. Passengers are differentiated by their assessed threat value, which is determined by a passenger prescreening system. Their objective is to maximize the total security, which is derived from the probability of detection in each class and the proportion of assessed threat value in the corresponding group of passengers. The MAP assigns passengers such that budget and passenger assignment constraints are not violated, but does not consider the number of screening devices to use for each passenger class. They also consider the case where the passengers are assumed to be identical (i.e., passenger prescreening is not used). The results show that, if the passenger prescreening system is highly effective in identifying passenger risk, the security system is more efficient when the passengers are distinguishable. They also found that two or three classes may be sufficient. As new technologies are expected to be available in the near future, updating security systems at the airport is a critical problem. There has been little research in the subject of allocating new screening devices for airport security systems that considers the probability of threat detection and device allocation costs. In this research, allocating both newly developed and existing devices used by an airport has been explicitly considered for an airport security system.

16

2.3 Complexity Theory Many optimization problems are NP-hard; hence, it is worthwhile to discuss complexity and NP-completeness theory. According to Hall (2000) “Complexity Theory” or “Computational complexity” has two meanings. One is a measure of how many steps are required for the algorithm to solve the problem in the worst case for an input of a given size. On the other hand, the theory of computational complexity provides a mathematical framework in which computational problems can be classified as “easy” or “hard” to solve. To begin, a problem is a general question requiring an answer, whereas an instance of a problem will include an exact value of the data. An algorithm is a finite sequence of steps guaranteed to find the correct solution to any instance. We can define, more mathematically, a computational problem to be a function h that maps an instance x in some given domain to an output h(x) in some given range. The complexity of an algorithm is an upper bound of f(n), the number of steps to solve the problem in the worst case given that the input size of the problem is equal to n. We shall say that f(n) = O(g(n)) if there is a positive number c and n0 such that, for every n ≥ n0 , f(n) ≤ cg(n). In addition, when f(n) = O(g(n)), we call g(n) an asymptotic upper bound for f(n). An algorithm is a polynomial-time algorithm (has a polynomial-time complexity) if its running time function f(n) = O(p(n)) for some polynomial function p(n). Problems with polynomial-time algorithms are said to be “easy” and those algorithms are considered efficient. For example, consider two algorithms with asymptotic upper bound equal to n2 (polynomial) and 2n (exponential), respectively. When input size n = 10, the first algorithm requires 102 = 100 steps in the worst case, while the second algorithm

17

requires 210 = 1024 steps. However, if the input size n increases to 100, the difference between the number of steps required for those two algorithms could be as large as 1030. Consequently, many mathematicians prefer problem that can be solved in polynomial time. The theory of NP-completeness is applied only to the decision problem class. This class consists of those problems having only two possible solutions, yes and no. We say that a decision problem belongs to class P if the solution can be obtained by a polynomial-time algorithm, where P stands for polynomial time. A particular superset of class P is class NP. The name NP stands for nondeterministic polynomial, and class NP refers to a set of problems with the following property: any “yes” instance can be verified in polynomial time. Examples of problems belonging to NP are Traveling Salesman, Subgraph Isomorphism, and Hamiltonian Cycle. A famous open problem about complexity is whether P is the same as NP or if P is a proper subset of NP. To study that problem, we introduce the most important subset of NP, called NP-complete. The problem П is said to be in NP-complete if П is in NP and any other NP problems are polynomial time reducible to П. We say problem A is polynomial time reducible to problem B, if there exists a polynomial time mapping function f with the following property: for every x, x ∈ A ? f(x) ∈ B . The theory of NPcompleteness was first introduced by Cook (1971). He states that if a polynomial time algorithm exists for any problem in NP-complete, then every problem in NP can be solved in polynomial time. In other words, P =NP if and only if one of the NP-complete problems is solvable in polynomial time. One more property of NP-completeness is that,

18

a decision problem Q1 is NP-complete if Q1 is polynomial time reducible to any other NP-complete problem. Gary and Johnson (1979) show that the following six problems are NP-complete and can be used as references: 3-Satisfiability, 3-dimensional matching, Vertex cover, Clique, Hamiltonian Circuit and Partition. Any optimization problem can be converted into an associated decision problem. For example, consider the maximization problem f(x), the associated decision problem is to answer the question “Is f(x) ≤ B?” for some number B. An optimization problem belongs to class NP-hard if its associated decision problem is NP-complete. The word hard is used to imply that any problem in NP-hard is at least as hard as any problem in NP.

2.4 Optimization Problems and Algorithms The objective of the model in this paper is to assign both the number and type of screening devices that optimize the probability of threat detection. The model shares certain characteristics with well-known optimization problems, such as the Integer Nonlinear Programming problem, (Nonlinear) Knapsack problem, Allocation problem, and Reliability problem. Many algorithms, heuristics and complexity issues have been derived for each type of problem. An analysis for each type depends on problem characteristics, such as the types of constraints and/or objective function, the types of decision variables, number of constraints, etc. For example, certain algorithms require objective functions and/or constraints to be concave, convex, differentiable, monotonic, separable, in quadratic form, etc. For some methods, the magnitude of the largest

19

numbers and the number of decision variables can exponentially increase computational time. We present a brief review of the Integer Nonlinear Programming problem, the Knapsack Problem, and the Reliability Problem in the following subsections.

2.4.1 Integer Nonlinear Programming Problems Since the objective function of the model in this research is nonlinear and the decision variables are integers, the model is an Integer Nonlinear Program (INLP). Complexity analysis shows that INLP-decision problems are NP-complete (Nemhauser and Wolsey, 1988). In general, INLPs are nonconvex, which implies the potential existence of multiple local solutions. For a problem with continuous variables, smooth, and nonconvex nonlinear programming, Murty and Kabadi (1987) show that verifying that a given feasible solution is not a local minimum and verifying the objective function is not bounded below on the set of feasible solutions are NP-complete problems. As a result, finding a global solution for an INLP problem is NP-hard. Despite the hardness of INLP problems, Floudas (1995) states that several algorithms and heuristics have been studied and developed to solve INLP problems, for examples, Generalized Benders Decomposition (GBD), Branch and Bound (BB), Outer Approximation (OA) , Generalized Outer Approximation (GOA), etc. Pardalos (1996) presents an overview of how to formulate certain integer (or combinatorial) optimization problems into equivalent continuous problems (relaxation problem) which can be used to obtain approximate algorithms or to compute sub-optimal solutions. Bretthauer and Shetty (1995) develop a Branch and Bound algorithm with a

20

general framework for solving a continuous relaxation of the problem. Heuristics and reoptimization procedures that round a continuous solution to obtain a feasible solution are also given to improve the performance of the Branch and Bound method. They also demonstrate how to modify the algorithm to handle nonconvex problems. Following their work in 1995, Bretthauer and Shetty (2002) present a “pegging algorithm” for solving a continuous relaxation of the problem and embed this algorithm in Branch and Bound. Computational results show that a pegging branch and bound algorithm is more efficient than a multiplier search branch and bound algorithm, dynamic programming, and a 0,1 linearization method. However, this method only applies to a differentiable convex objective function over a single differentiable convex constraint and bounded integer variables. Lawler and Belle (1966) propose a simple and effective search method for solving discrete optimization problems with a monotonic objective function and arbitrary constraints. Basically, for binary optimization problems, the method arranges all possible (not necessarily feasible) solutions from lowest to highest according to their numerical values, given by n( x) = x1 2 n ?1 + x2 2 n ? 2 + ... + xn 2 0 . Instead of searching all possible solutions, they provide three rules to skip searching many possible solutions, thus reducing the enumeration steps. Nevertheless, Lawler and Belle’s method cannot handle problems with quadratic objective functions or non-separable objective functions. Moreover, Lawler and Belle’s algorithm requires all decision variables to be binary. Integer variables must be transformed into binary variables before applying an algorithm, thus increasing the size of the system and computational time. Mao and Wallingford

21

(1968) extend Lawler and Belle’s method to handle problems with a quadratic objective functions by modifying the three skipping rules in the enumerative algorithm.

2.4.2 Knapsack Problems A Knapsack Problem is a special type of integer program. In general, the problem is to select a collection of items that do not exceed the capacity of a knapsack such that the total value is as large as possible given the weight and value of each item. The Knapsack problem is also NP-complete (Lueker, 1975). Kellerer et al. (2004) show that the following special cases of decision knapsack problems are NP-complete: the bounded knapsack, the unbounded knapsack, the unbounded subset sum, the multiple knapsack, the multidimensional knapsack, the multiple-choice, and the quadratic knapsack problems. Even though the Integer Knapsack Problem is NP-complete, Lusena (1999) shows that the time to solve the problem is polynomial in the size of the instance plus the capacity of the knapsack using a dynamic programming approach. In addition, the integer knapsack problem is solvable in polynomial time in the size of the instance plus the sizes of the items. Marthur et al. (1986) suggest solving the nonlinear knapsack problem (NLK) by using a piecewise linear approximation of a concave and increasing objective function, and convex and decreasing constraints, to convert NLK into a classical knapsack problem. As a result of the approximation, an efficient search enumeration procedure can be developed by tightening the bounds of the decision variables. Kovalyov (1996)

22

presents an approximation algorithm for NLK, without restriction to concave or convex objective and constraint functions. His algorithm is claimed to solve the problem in O nU ∑ j =1 min{nU /(εL), b j } /(εL) , where L and U are a lower and upper bounds for the

n

(

)

optimal solution, respectively.

2.4.3 Reliability Problems Detecting threats in airport security systems relies on a configuration of explosive screening devices. The concept of reliability optimization problems can be applied to maximizing the probability of threat detection. There are many types of configurations when considering system reliability, for examples, Series, Parallel, Series-Parallel, Standby, Bridge network, and Nonseries-nonparallel. In this research, the probability of threat detection for a combination of explosive screening devices is considered to be a parallel system, because the threat detection system fails when all of the explosive screening devices fail. The reliability optimization problem can be formulated as an integer programming problem, since the number of components must be integer valued to make any sense. Many approximate and exact techniques have been applied to optimize the reliability of a system. It is obvious that making use of component redundancy increases system reliability. For large scale reliability problems, an approximate approach is preferable to an exact algorithm, since it is faster and less costly. Sharma and Venkateswaran (1971) propose an intuitive heuristic approach for allocating redundancy among subsystems. The procedure is, basically, to add one component to the stage with

23

the lowest reliability. Nevertheless, the approach may not yield an optimal solution if the stages have components of similar reliability but extremely different cost. Revising the work of Sharma and Venkateswaran (1971), Aggarwal et al. (1975) make use of slack variables to prevent choosing only one component from the lowest reliability stage. This approach adds a component based on a ratio of “increment increases in reliability” to the “product of increment increases in resources usage.” Misra (1971) claims that the method of Lawler and Belle (1966) is easy to implement and economical to solve many types of reliability problems. Examples of the system in his paper include problems that have linear objective functions, nonlinear objective functions and linear constraints, and nonlinear objective functions and nonlinear constraints. Misra and Sharma (1973) develop a zero-one programming algorithm to solve a series of reliability systems. Their technique depends on the nonbinary tree-search, where the upper bound is obtained using graph theory. However, the reliability system problem must have at least one component in each stage (or subsystem). Sharma and Misra (1990) propose another search method which can solve problems independent of the nature of the objective function. That is, it can be applied to problems with following properties: series, parallel or mixed, and non-separable, integer and/or binary decision problems (as opposed to Lawler and Belle’s algorithm). Moreover, the problem objective function is not required to have a partial derivative with respect to the components.

24

In the next chapter, we formulate the explosive scanning device allocation problem for a single airport station model based on reliability and integer programming problems. The model is then proved to be NP-Hard by showing a polynomial transformation from the Binary Knapsack problem. In Chapter 4, we solve the single station model by applying Lawler and Bell’s partial enumeration algorithm (1966).

25

CHAPTER III EXPLOSIVE SCREENING DEVICE ALLOCATION MODEL FOR A SINGLE AIRPORT

3.1 Introduction In this chapter, the airport security system as well as the screening technology allocation process are described. The objective of this research is to develop a screening device allocation model for airport security systems. The goal is to maximize the expected total security. There are three primary means passengers can use to get a threat onto an aircraft; on their person, in their carry-on baggage or in their checked baggage. Carry-on baggage refers to items that a passenger carries onto the plane other than on their person. It includes luggage, briefcases, purses, computer bags, coats, etc. During the check-in process, a passenger prescreening system classifies passengers into different classes (for example, selectee and non-selectee classes) according to the level of threat they pose to an aircraft. Passengers, checked baggage, and carry-on baggage are screened by a predefined combination of screening devices or procedures, which search for threats. Each screening device has different characteristics, such a processing rate, operating and installation costs, and the probability of detecting a threat given that there is a threat. If many types of detection devices are required for a particular class of passengers, the security system can be costly, and passengers may experience an excessive delay. On the other hand, too few detection devices can cause a great danger to

26

an aircraft as the system might not be able to detect threats from a terrorist. The challenge is to assign appropriate types of screening devices to each class and to determine the number of the devices needed in order to maximize the total security as well as ensure that the cost constraint is satisfied. Different types of technologies are typically used to screen the different means of threat; passengers, carry-on baggage, checked baggage. If a specific technology is used, several devices are needed because items are screened at different physical locations. As a result, the problem of allocating sets of technologies and devices for the different means can be separated into three independent sub-problems. This research focuses on the problem of grouping technologies for screening carry-on baggage. The resulting approach can also be applied to the problems of allocating technologies to groups of checked baggage or groups of passengers.

3.2 Single Airport Scenario In practice, passengers are classified into different classes by the passenger prescreening system. It is assumed that the passenger prescreening system is highly effective in distinguishing passengers, i.e. high-risk passenger receives high assessed threat value, and low-risk passenger receives low assessed threat value. The carry-on baggage within theses classes is then inspected by different combinations of screening devices and/or procedures. If one of the screening devices or procedures detects a threat, the bag receives extra scrutiny. The bag is allowed onto an aircraft only if it is cleared by all of the assigned screening devices. The expected number of bags as well as an average

27

assessed threat value are given for each class of passengers. The costs associated with the screening devices include fixed, operating and installation costs. The needed screening capacity can be obtained by using existing devices or purchasing new devices. The installation cost occurs only if new devices are installed. There is no installation cost for using existing devices. The model also includes the situation where there is new technology available. T can determine whether to deploy a new technology should be deployed (probably faster and more accurate but more expense to purchase and install) or existing device should be retained (less costly but maybe slower and less accurate). Each device has different probability of threat detection. This research assumes that each device is independent, i.e. installing one device does not affect the performance and cost of other devices. For example, the probability of threat detection of device type j given that it is a threat, Pj, is equal to the probability of threat detection of device type j given that it is a threat and the threat is not detected by the prior device, or, Pj, = P{Device j detects a threat | threat is present} = P{Device j detects a threat | threat is present & device j-1 cannot detect a threat} = P{Device j detects a threat | threat is present & device j-1 is not employed}.

Generally, baggage will receive a secondary screening if it sets off the alarm of any screening devices. In this research, the focus is on the primary screening procedure and the model is assumed that there are enough resources for resolving the alarm by the secondary procedure. Therefore, this research assumes that there is no false alarm in the system, i.e. when the primary device detects a threat, it is actually a threat.

28

Several assumptions are made regarding the single airport scenario. They are summarized as follows: i. A passenger prescreening system is used to assign passengers (and their baggage) to classes. ii. The assessed threat value (assigned by a passenger prescreening system) for each passenger is the accurate representation of that passenger’s true level of risk. iii. iv. v. vi. Each device is independent. The quantity of new devices is unlimited. There is no cost associated with uninstalling existing devices. There is no false alarm.

The objective of the model is to assign both types and numbers of screening devices and/or procedures to each group so that the total probability of detection is maximized subject to cost and capacity constraints. We consider assigning a combination of screening devices to maximize the security system given classes f passengers assigning. Compare to previous work, Rausch (2002) does not take the probability of threat detection into account, only device capacity and associated costs. Lazar Babu (2003) randomly assigns passengers to classess, so passengers are indistinguishable, which may not be realistic given passenger prescreening information. McLay et al. (2006) focus on grouping passengers, not forming combinations of devices, and

29

screening device capacity is not considered. Table 3.1 outlines the characteristic that distinguish our model from those studies in the past.

Table 3.1. Model Comparison. Single Airport Model Maximize total security Number and Type of screening device for each class M classes of passengers with average assessed threat value Yes Rausch (2002) Maximize throughput Number of device for each airport McLay (2006) Maximize total security Lazar Babu (2003) Minimize false alarm Passenger grouping and number and type of screening device for each group Number of passengers

Objective

Decision

Passenger Grouping

Predefined

Production schedule

Assessed threat value of each passenger and combinations of devices Yes

Passengers are distinguishable Consider probability of threat detection

No

No

Yes

No

Yes

Yes

30

3.3 Notation This section presents the notation used in the single airport model. M – number of passenger classes, D – number of screening devices types, Ai – average assessed threat value of passengers in class i (i = 1, 2, …., M), Bi – the number of bags in class i (i = 1, 2, …., M), BT – total number of bags =

∑B

i =1

M

i

,

Fj – fixed cost ($/device) of device type j (j = 1,2, …., D), Ij – installation cost ($/device) of device type j (j = 1,2, …., D), Oj – operating cost ($/bag) of device type j (j = 1, 2, …., D), Cj – capacity (bags/hour) of device type j (j = 1, 2, …., D), Ej – a number of existing devices of device type j (j = 1, 2, …., D), Pj – probability of detecting a threat given that there is a threat for device type j (j = 1, 2, …., D), TB – total budget ($).

Decision Variables: Xij – binary variable where Xij =1 if device type j is used to screen class i bags, Xij = 0 otherwise, Yj – number of devices of type j required (integer).

31

3.4 Single Airport Model The allocation model needs to assign both numbers and types of devices to be used for each class. Once the number of a particular device type required is found, the number of devices to be installed must be determined. We can determine the number of device type j to be installed as the number of device type j required minus the number of device type j currently existing at the airport. Denote the number of device type j to be installed by Sj , hence Sj = Yj – Ej, for j = 1, 2,…, D. Since the uninstallation cost is not considered here, Sj = Yj – Ej if Yj > Ej, and Sj = 0 otherwise, for j = 1, 2,…, D. Subsequently, the total installation cost can be calculated as (3.1)

∑I S

j =1 j

D

j

and the total fixed cost can be calculated as

∑F Y

j =1 j

D

j

,

where Ij and Fj are estimated hourly installation and fixed costs, respectively, and are independent of the volume of checked bags screened. Fixed cost includes annual maintenance and annual lease expenses. Installation cost includes purchasing and preparation costs. These costs are obtained based on the expected useful life of the device and the number of operational hours per year (per day times 365).

32

An operating cost for device type j depends on the number of bags in each class screened by device type j. It is calculated as

∑O B X

i =1 j i

M

ij

.

Therefore the total operating cost is given by

∑∑ O B X

j =1 i =1 j i

D

M

ij

,

where Oj is expected cost of operating a device type j, based on the salaries of the screening device operators. For a total budget, TB, the cost constraint is

∑ (F Y

D j =1 j

j

+ I j S j ) + ∑∑ O j Bi X ij ≤ TB ,

D M j =1 i =1

(3.2)

where the first term consists of fixed and installation costs, and the second term represents the operating cost. It is desirable to determine the number of screening devices required to handle the number of bags. Therefore the capacity of screening devices should be considered. For each device type j, the total number of bags screened is

∑B X

i =1 i

M

ij

.

Then, the total number of device type j needed is the total number of bags screened by device type j divided by capacity of device type j,

∑B X

i =1 i

M

ij

Cj .

33

However, the number of devices is an integer. The above equation will not necessarily yield integer values, so we set the lower bound by defining the following capacity constraints:

Yj ≥

∑B X

i =1 i

M

ij

Cj

M

,

or

C jY j ≥ ∑ Bi X ij for j = 1, 2,…, D.

i =1

(3.3)

Similarly the upper bound capacity constraints are:

C j Y j < ∑ Bi X ij + C j , for j = 1, 2,…, D.

i =1 M

(3.4)

On the left-hand side, C jY j is the total capacity of all type j devices, the right-hand side is the total requirement for device j (sum of the numbers of bags class i which are required screened by device j). The lower bound capacity constraints (3.3) ensure that the number of device j is sufficient to handle number of bags to be screened by device j. The upper bound capacity constraints (3.4) guarantee that when device k is not used for any class of passenger, i.e. Xik = 0 for all i, Yk is also equal to zero. The probability of threat detection for each class is defined as the probability that one of the device types in the combination can detect a threat. To simplify, denote Pj as the probability that a threat is detected by device type j given that a threat is present, and let Li denote the probability that a threat is detected by a combination of devices for class i given a threat is present, then it follows from the independence assumption that

34

Li =

P{a threat is detected by one of the devices screening class i passengers | a threat is present}

= 1-P{threat is not detected by any of the devices in class i passengers | a threat is present} Li = 1 ? ∏ (1 ? X ij Pj ), for i = 1, 2, …., M.

D j =1

(3.5)

An example is illustrated in Figure 3.1, with three explosive screening device types for class i. A bag is required to pass through all three screening device types before being carried on to an aircraft. A bag is taken out if one of the device types detects a threat.

Threat detected (Xi3P3)(1- Xi1P1)(1- Xi2P2) (Xi1P1) (Xi2P2)(1- Xi1P1) Boarding (threat not detected) Device1 (1- Xi1P1 ) Device2 (1- Xi1P1)(1- Xi2P2) Device3 (1- Xi1P1)(1- Xi2P2)(1-Xi3P3)

Passenger with a threat

Figure 3.1. Example of a Combination of Three Device Types.

35

We assume that the number of bags, Bi, and the average assessed threat value of passengers in each class, Ai, are known. The risk level of class i is defined as the proportion of the total assessed threat value assigned to class i,

Ri = Ai Bi

1

∑AB

i =1 i

M

.

(3.6)

i

Alternatively, risk level is the conditional probability that class i contains a threat given that the system contains a threat. To be precise, define the following incidents: Ti = class i contains at least one threat, ST = the system contains at least one threat. Then the risk level of class i, Ri, can also be represented as P(Ti | ST). The true alarm rate of a single screening device is the probability that the device detects a threat given that the bag contains a threat. Similarly, the true alarm rate of a combination of device types for class i is the probability that at least one of the device types can detect a threat given that class i contains a threat. Subsequently, the total security is defined as the probability that a threat is detected given that there is a threat in the system. To determine the total security, let Di = a threat is detected in class i passengers (i.e. a class i bag contains a threat and the threat is detected given that the system contains a threat). Therefore, P(Di ) can be calculated as (the probability that class i contains a threat) x (the probability a threat is detected by a combination of device types for class i) or P(Di ) = P(Ti | ST) Li = Ri Li. (3.7)

36

As a result, the total security is

∑ P( Di ) =∑ Li Ri .

i =1 i =1

M

M

(3.8)

The allocation model can be formulated as an integer program, with binary decision variables Xij = 1(0), if device type j is (not) used in class i, and integer decision variables Yj = the number of screening device type j required, Yj ≥ 0. Maximize subject to

D

∑L R

i =1 i

j j G

M

i

∑ (Y F

j =1

+ S j I j ) + ∑∑ X ij O j Bi ≤ TB

M D i =1 j =1

Y j C j ≥ ∑ (X ij Bi ) , j = 1,2, . . ., D Y j C j < ∑ Bi X ij + C j , j = 1,2, . . ., D Xij ∈ {0,1} Yj ∈ Ζ +

1 and Li = 1 ? ∏ (1 ? X ij Pj ).

D j =1

i =1 i =1 M

where

Ri = Ai Bi

∑AB

i =1 i

M

i

The objective function maximizes the total security. The first constraint is the budget constraint. The second constraint ensures that the screening device has sufficient capacity to handle the number of bags. The third constraint ensures capacity is not assigned if not needed. The last three constraints restrict integer and binary decision variables.

37

3.5 Computational Complexity of the Model Many optimization problems are proven to be NP-hard. Practically, this means the number of operations required to find the optimal solution in the worst case is exponential (not polynomial) in the input length. Finding the optimal solution for a problem with a large input size could take centuries even with the fastest computer. Since, it is unlikely that there is an exact algorithm for solving NP-hard problems in polynomial time, we must focus on finding an acceptable solution instead of an optimal solution. In this section, we prove that the explosive detection device allocation model is an NP-hard problem. We can consider a problem П to be NP-hard, if any problem in NP can be transformed to П in polynomial time, whether or not П is a member of NP. One approach to proving that a problem П is NP-hard is to show that there exists an NPcomplete problem П? that is polynomial time reducible to П, i.e. if П?∈NP-complete and П? ≤p П, then П is NP-hard, where П? ≤p П means П? is polynomial time reducible to П. Another technique proposed by Gary and Johnson (1979) is to restrict a problem П to be identical to a known NP-complete problem П?. In other words, we have to show that П contains a П? as a special instance of the problem. The Binary Knapsack Problem or 0-1 Integer Programming Feasibility Problem is known to be NP-complete problem (Nemhauser and Wolsey 1999). The Binary Knapsack Problem (BKP) is

INSTANCE:

? ?

m∈Z+ — number of items, uj ∈{0,1} — binary variables, j = 1, 2, …., m,

38

? ? ? ?

m

vj ∈R+ — a value of each item, wj ∈ R+ — a weight of each item, b∈ R+ — a weight constraint, k∈ R+ — a value goal.

m

QUESTION: Is there a u = (u1, u2, u3,..., um) such that

∑ (u w ) ≤ b and ∑ (u v ) ≥ k , where uj ∈ {0,1}, j = 1, 2,…, m?

j =1 j j j =1 j j

(3.9)

The decision problem for the explosive screening device allocation model can similarly be defined as

INSTANCE:

? ? ? ? ? ? ? ? ? ? ? ? ?

D — number of device types, M — number of classes, Bi ∈Z+ — a number of bags for class i, Ai ∈ R+ — an average assessed threat value for each class i passenger, Ej ∈Z+ — a number of existing machine for each device type, Cj ∈ R+ — a capacity of device type j, β ∈ [1, ∞] — an upper bound coefficient, Fj ∈ R+ — a fixed cost of device type j, Ij ∈R+ — an installation cost of device type j, Oj ∈R+ — an operating cost of device type j, Pj ∈(0, 1] — the probability of threat detection given that there is a threat for each device type, K∈R+ — a value goal, TB∈R+ — a total budget.

? X 11 ?X ? 21 ? : ?X ? M1 ? X 12 X 22 : XM2 ... ... : ... X 1D ? X 2D ? ? : ? X MD ? ? ?

QUESTION:

Is there X =

and Y = (Y1, Y2, Y3,..., YM) such that

D ? ? ∑ ?1 ? ∏ (1 ? X ij Pj )?Ri ≥ K , i =1 ? j ? M

(3.10) (3.11)

∑ (Y F

D j =1 j

j

+ S j I j ) + ∑∑ (X ij O j Bi ) ≤ TB ,

M D i =1 j =1

39

Y j C j ≥ ∑ (X ij Bi ), for j = 1,2, . . . , D,

M

(3.12) (3.13)

Y j C j < ∑ Bi X ij + C j , for j = 1,2, . . . D , Xij ∈ {0,1}, Yj ∈ Ζ + , where Ri = Ai Bi 1 for i = 1, 2, …., M ,

i =1

i =1 M

∑AB

i =1 i

M

i

and

Sj = Yj – Ej if Yj > Ej, and Sj = 0 otherwise for j = 1, 2,…, D?

Theorem 3.1 The decision problem for the explosive screening device allocation model is NP-hard. Proof: The proof consists of two parts. Part I shows the polynomial transformation from the BKP decision problem to the explosive scanning device allocation decision problem, and Part II shows one-to-one correspondence between the two problems. (Construction suggested by Dr. Christopher J. Monico.) Part I: Transformation From an arbitrary instance of the Binary Knapsack Decision Problem we can construct a particular instance of the Explosive Screening Device Allocation Decision Problem by a polynomial transformation as follows: i. D = m, number of device types equals number of items of BKP.

40

ii.

M = 1 (i.e. only one group of passengers), then Bi = 0 and Ai = 0 for i ≥ 2, and since (3.6) Ri = Ai Bi 1 , we have R1= 1.

∑AB

i =1 i

M

i

iii.

Ej = 0, then from (3.1) Yj = Sj for j = 1, 2,…, D. (i.e. there is no existing device available).

iv.

Cj = B1 for j = 1, 2,…, m (i.e. the capacity of each machine is equal to the number of bags in class 1).

v. vi.

β = 1, then from (iv) and capacity constraint (3.12) and (3.13), Yj = X1j. from equation (3.11),

∑ (Y F

D j =1 j

j

+ Y j I j ) + ∑∑ (X ij O j Bi ) , since we have M = 1

M D i =1 j =1

and D = m, (3.11) becomes

∑ (Y F

m j =1 j

j

+ Y j I j ) + ∑∑ (X ij O j B j ) = ∑ (X 1 j (O j B j + I j + F j ) ) .

1 m m i =1 j =1 j =1

?v j

vii.

viii.

Pj = 1 ? e

, K = 1 ? e ? k , TB = b, and F j + I j + O j B1 = w j .

Then from (i) to (vii), we have the problem in the form of the explosive screening device allocation decision problem as

m ? ? ?v 1 ? ∏ 1 ? (1 ? e j ) X 1 j ? ≥ 1 ? e ? k , ? j =1 ? ?

(

)

∑ (X

m j =1

1j

w j ) ≤ b and Xij ∈ {0,1}.

(3.14)

The transformation requires O(m) time, and note that 1 ? e value since 0 ≤ 1 ? e

?v j

?v j

is a valid probability

≤ 1.

41

Part II one-to-one correspondence ix. To show that the “yes” (or “no”) solution of the Binary Knapsack problem is also a “yes” (or “no”) solution of the decision problem for explosive screening device allocation model, suppose that a binary vector u = (u1, u2, u3,..., um) is a “yes” solution of the BKP, i.e. u satisfies

∑ (u w ) ≤ b

m j =1 j j

(3.15)

and x.

∑ (u v ) ≥ k

m j =1 j j

(3.16)

From (vii), w j = F j + I j + O j B1 , therefore (3.15) becomes

∑ (u

m j =1

j

( F j + I j + O j B1 ) ) ≤ b . Since b = TB, we have ∑ (u j ( F j + I j + O j B1 ) ) ≤ TB .

m j =1

Since D = m, then u also satisfies the budget constraint of the explosive screening device allocation decision problem. xi. To show that u that satisfies objective constraint (3.15) also satisfies objective constraint of the explosive scanning device allocation problem (3.10), from (vii),

Pj = 1 ? e

?v j

, then vj = (? ln(1 ? Pj ) ) , therefore (3.16) becomes

∑ (u (? ln(1 ? P ))) ≥ k .

m j =1 j j

Since ui is a binary variable, u j (? ln(1 ? Pj ) ) = (? ln(1 ? Pj u j ) ) , then we have

∑ (? ln(1 ? P u )) ≥ k ,

m j =1 j j

42

? ∑ (ln(1 ? Pj u j ) ) ≥ k ,

m j =1

∑ (ln(1 ? P u )) ≤ ?k .

m j =1 j j

Because the exponential function is monotone increasing, we have

? m ? exp? ∑ (ln(1 ? Pj u j ) )? ≤ exp(? k ) , ? ? ? j =1 ?

∏ exp(ln(1 ? P u )) ≤ exp(?k ) ,

m j j j =1

∏ ((1 ? P u )) ≤ exp(?k ) ,

m j j j =1

? ∏ ((1 ? Pj u j ) ) ≥ ? exp(? k ) ,

m j =1

1 ? ∏ ((1 ? Pj u j ) ) ≥ 1 ? exp(?k ) .

m j =1

Since K = 1 ? e ? k and D = m, u also satisfies the objective constraint of the explosive scanning device allocation problem. It follows that u is also the “yes” solution of the explosive screening device allocation decision problem. xii. To show that the “yes” (or “no”) solution of the decision problem of the explosive screening device allocation model is also a “yes” (or “no”) solution of the Binary Knapsack problem, suppose that a binary vector x = (x11, x12, x13,..., x1m) is a “yes” solution of the explosive screening device allocation decision problem, i.e.

x satisfies

43

[ 1 ? ∏ (1 ? x1 j Pj ) ] ≥ K

D j =1

(3.17)

and xiii.

∑ (x

D j =1

1j

( F j + I j + O j B1 ) ) ≤ TB .

(3.18)

From (vii) F j + I j + O j B1 = w j , then (3.18) becomes

∑ (x

D j =1

1j

w j ) ≤ TB = b

since m = D, then x also satisfies the weight constraint of BKP. xiv. Without loss of generality, assume that 0 < Pj < 1, then it is obvious that

K ∈ (0,1) . Therefore, we can transform (3.17) as follows,

D ? ? ?1 ? ∏ (1 ? x1 j Pj )? ≥ K , j =1 ? ?

∏ (1 ? x P ) ≤ 1 ? K .

D 1j j j =1

Since the logarithmic function is strictly increasing and both sides of the equation are greater than zero,

? D ? ln ? ∏ (1 ? x j Pj )? ≤ ln(1 – K), ? ? ? j =1 ?

∑ ln(1 ? x P ) ≤ ln(1 – K).

D j =1 1j j

(3.19)

Because xij is a binary variable, it is obvious that ln (1 ? x1 j Pj ) is equal to x1 j ln (1 ? Pj ) . Then (3.19) becomes

44

∑x

j =1

D

1j

ln (1 ? Pj ) ≤ ln(1 – K)

or equivalently

∑ x (? ln(1 ? P )) ≥ -ln(1 – K).

D j =1 1j j

(3.20)

From (vii), substituting Pj = 1 ? e

?v j

and K = 1 ? e ? k in (3.20), yields

?v j ?k

∑x

j =1

D

1j

(? ln(1 ? 1 + e )) ≥ ? ln(1 ? 1 + e

∑x

j =1 D 1j

),

(? ln(e )) ≥ ? ln(e

?v j

?k

),

∑x

j =1

D

1j

vj ≥ k .

Thus x also satisfies the objective constraint of BKP. xv. The one-to-one correspondence between the instances that preserve “yes” and “no” solutions is proven to be valid. That is, the “yes” (or “no”) solution of the Binary Knapsack problem is also the “yes” (or “no”) solution of the decision problem of explosive screening device allocation model and vice versa. Therefore it can be concluded that the explosive screening device allocation problem is NP-hard. ■

The explosive screening device allocation problem is formulated as an integer nonlinear optimization problem with two types of decision variables, integer and binary. The objective of the model is to assign a combination of screening devices that maximizes the probability of detecting a threat. The total budget and throughput of the system are taken into account. The model is proved to be NP-Hard by constructing a 45

polynomial transformation from a known NP-Complete problem; the Binary Knapsack problem. Therefore, a large scale explosive scanning device allocation model may not be solvable to optimality within a reasonable computational time. However, for a small size of the problem, i.e. three classes of passengers and seven devices, the model can be solved optimally in a reasonable time by a method of partial enumeration, which is described in Chapter 4.

46

CHAPTER V MULTIPLE AIRPORT PROBLEM

5.1 Introduction In Chapter 3, we presented the development of the optimization model for a single airport (station) security system, where the focus is on carry-on baggage screening. The model goal is to maximize the probability of threat detection while considering accounts for capacity and budget constraints. In Chapter 4, we discussed Lawler and Belle’s partial enumeration algorithm (Lawler and Belle, 1966) which is used to solve the single-station problem. The model is modified to be consistent with Lawler and Belle’s standard form. From the result of the experiments, it can be seen that Lawler and Belle’s algorithm is capable of solving a single station problem within a reasonable time. In order to prevent a threat from reaching an aircraft, we must consider a set of important commercial airports instead of just one airport, because terrorists can attempt to fly from one of many commercial airports. Moreover, new technology for detecting a threat is likely to be available in the near future, as the terrorists never stop looking for a way to breach the current security system. To prevent the newly emerging threats, the security system of all important airports must be renovated from time to time. When a new screening system becomes available, it will only be available in limited quantities; one of the most critical decisions is to determine which airports will receive the system first, and what to do with the existing systems. To assign the type and the number of scanning devices to a set of considered airports is very challenging. In order to overcome

66

such problems, in this chapter, the single airport model is expanded to include multiple airports in one optimization problem. The resulting model can be applied to a real-world problem when allocation of a new set of devices to a set of considered airports is required.

5.2 Multiple Airport Scenario In this chapter, we extend the number of airports considered in the optimization model. In the Multiple-Airports scenario, each airport has particular characteristics, such as the number of passengers, number of passenger classes, number and type of existing devices, and budget availability. Some airport assumptions are as presented in Chapter 3, and are summarized below. For each airport, ? ? Passengers are grouped according to their posed level of threat. Each group of passengers has to be inspected by a different combination of the screening devices or procedures. ? If one of the screening devices or procedures detects a threat, the bag receives extra scrutiny. The bag is allowed onto an aircraft only if it is cleared by all of the assigned screening devices. ? The expected number of bags as well as an average assessed threat value is given for each class of passengers. ? The needed screening capacity can be obtained by using existing devices or installing new devices.

67

?

There is no installation cost for using existing devices. In the multiple airports scenario, the costs associated with the screening devices,

including fixed, operating and installation costs are assumed to be the same for all airports. The difference between the single airport model and the multiple airports model is there will be limited resources for new devices and these must be allotted among the airports under consideration. Revised assumptions for the multiple airports problem are summarized as follows: i. A passenger prescreening system is used to assign passengers (and their baggage) to classes. ii. The assessed threat value (assigned by a passenger prescreening system) for each passenger is the accurate representation of that passenger’s true level of risk. iii. Each device is independent (i.e. installing one device does not affect the performance and cost of other devices). iv. v. vi. vii. The quantity of new devices is limited. The costs associated with the screening devices are the same for all airports. There is no cost associated with uninstalling existing devices. There is no false alarm.

5.3 Multiple Airport Model In this section, the development of the multiple airports model is presented. We start by defining the new notation used in the multiple airports model in Section 5.3.1,

68

and creating the new objective function and additional constraints for the multiple airports optimization problem in Section 5.3.2.

5.3.1 Notation For T airports considered, we add another index, k, to represent the kth airport, where k = 1, 2, . . ., T. We redefine the notation used in Chapter 3 as, Mk – number of passenger classes at airport k, Aik – average assessed threat value of passengers in class i at airport k, Bik – the number of bags in class i at airport k, D – number of screening devices types, Fj – fixed cost ($/device) of device type j, Ij – installation cost ($/device) of device type j, Oj – operating cost ($/bag) of device type j, Cj – capacity (bags/hour) of device type j, Eik – the number of existing devices of device type j at airport k, Pj – probability of detecting a threat given that there is a threat for device type j, TBk – total budget($) for airport k, Uj – Number of device type j available. Decision Variables: Xijk – binary variable where Xij =1 if device type j is used to screen class i bags at airport k, Xijk = 0 otherwise, Yjk – number of device type j required at airport k (integer).

69

5.3.2 The Model From Equation (3.5) of the single airport model, for class i passengers at airport k, the probability of threat detection of a combination of devices given that there is a threat is Lik = 1 ? ∏ (1 ? X ijk Pj ) .

D j =1

(5.1)

We assume that the assessed threat value for each passenger is the accurate representation of that passenger’s true level of risk, and define Rik as the proportion of the total assessed threat value assigned to class i of airport k or,

Rik = Aik Bik

1

∑A

i =1

Mk

. Bik

(5.2)

ik

Therefore, a measure of security for airport k is

∑L

i =1

Mk

ik

Rik .

(5.3)

However, with this measure, each airport will have the same weight, and we will not be able to distinguish the important airports from other airports. We can develop a measure of performance for the multiple airports model by weighting each airport by the number of bags to be screened at that airport. Then, the new measure for each airport, namely security level (SLk), is a sum of the probability of threat detection times the number of bags in that class

SLk = ∑ Bik Lik Ri k ,

i =1

Mk

(5.4)

70

and the total security level of all airports is

∑ SL

k =1

T

k

.

(5.5)

The total cost of airport k can be calculated in the same manner as the single airport model,

∑ (Y

D j =1

jk F j + S jk I j ) + ∑∑ X ijk O j Bik , Mk D i =1 j =1

where the number of device type j to be installed, Sjk = Yjk – Ejk if Yjk > Ejk, and Sjk = 0 otherwise, for j = 1, 2,…, D. Then, the cost and capacity constraints for airport k can be similarly defined as

∑ (Y

D j =1

jk F j + S jk I j ) + ∑∑ X ijk O j Bik ≤ TBk , Mk D i =1 j =1

(5.6)

Y jk C j ≥ ∑ (X ijk Bik ) , j = 1,2, . . ., D,

Mk i =1

(5.7)

Y jk C j < ∑ (X ijk Bik ) + C j , j = 1,2, . . ., D.

M i =1

(5.8)

As the availability of each device type is limited, we need to add one constraint for each type of the device, for j = 1,2, . . ., D,

∑S

k =1

T

jk

≤U j.

(5.9)

71

The explosive scanning device allocation problem for multiple airports problem is formulated as Maximize

∑ SL

k =1

T

k

(5.10)

M D

subject to

∑ (Y

D j =1

jk F j + S jk I j ) + ∑∑ X ijk O j Bik ≤ TBk , for k = 1,2, . . ., T, i =1 j =1

G

(5.11) (5.12) (5.13) (5.14)

Y jk C j ≥ ∑ (X ijk Bik ) , for j = 1,2, . . ., D, and k = 1,2, . . ., T, Y jk C j < ∑ (X ijk Bik ) + C j , for j = 1,2, . . ., D, and k = 1,2, . . ., T,

i =1 i =1 M

∑S

k =1

T

jk

≤U j

for j = 1,2, . . ., D,

Xijk ∈ {0,1}, Yj ∈ Ζ + , where SLk = ∑ Bik Lik Ri k ,

i =1 M

? ? Rik = ? Aik Bik ? ? ?

? ? D 1 ? and Lik = 1 ? ∏ (1 ? X ijk Pj ). Mk ? j =1 ∑ Aik Bik ? i =1 ?k

The objective function is the sum of security levels which is the product of number of bags in the class, probability of threat detection, and the proportion of the total assessed threat value assigned to that class. Constraint (5.11) limits the cost of each airport to be within that airport’s budget. Constraint (5.12) ensures that the screening device has enough capacity to handle number of bags in that airport. Constraint (5.13)

72

guarantees that the device is not assigned if it is not needed. Constraint (5.14) limits the total number of new devices to be less than or equal to the number of new devices available. The problem becomes a large scale optimization problem as the number of airports increases. It can be seen that, the overall objective function which is to be maximized is the sum of each airport’s performance function, i.e. the objective function is separable. It is obvious that this problem has a block structure: constraints (5.11), (5.12) and (5.13) are related to separate airports, which can be considered as block constraints for each airport (Tsurkov, 2001). However, there are interactions among separate airports (or blocks) in constraint (5.14). In other words, constraint (5.14) connects separate blocks in the model, and these constrains will be referred as binding constraints of the multiple airports model for the rest of this dissertation. With the existence of the binding constraints, it is clear that the model cannot be decomposed into independent single-airport problems which can be independently optimized. The problem also belongs to the NP-hard class, since the restricted version, T = 1 (one airport) was shown to be NP-hard problem in Section 3.5. Therefore, for a large number of airports, the problem may not be solvable to optimality in a reasonable time. We must focus on finding an acceptable solution by using a heuristic instead of finding an optimal solution by an exact algorithm. Proposed methodologies are discussed in the next chapter.

73

CHAPTER IV METHODOLOGY AND EXAMPLES FOR SINGLE AIRPORT MODEL

4.1 Introduction The explosive screening device allocation model can be considered as a Nonlinear Integer Problem (NIP), because the objective function, maximize total security, is nonlinear and the decision variables, the number of screening devices and use/not use binary variables are integers. NIP usually belongs to the class of NP-Complete problems, which means that finding an optimal solution for a large problem may not be possible in a reasonable time. Many techniques have been proposed to solve Nonlinear Integer Programming problems, for example, Greedy Heuristics, Dynamic Programming, Branch-and-Bound, Outer Geometric Programming, Method of Lagrange Multipliers, etc. Unfortunately, there is no unique algorithm that efficiently solves every Nonlinear Integer Programming problem. The performance of each method depends on problem characteristics, such as the types of constraints and/or objective function, the types of decision variables, number of constraints, etc. For example, certain methods require the objective function to be concave, convex, differentiable, in quadratic form, etc. For some methods, computational time depends on the magnitude of the largest parameter. In this chapter, a method for solving the explosive screening device allocation problem is presented using the Partial Enumeration Method proposed by Lawler and

47

Belle (1966). A brief review of the Lawler and Belle (LB) method is provided in Section 4.2, followed by a modification of the explosive screening device allocation model in Section 4.3, lastly, an example of the model and results are presented in Section 4.4.

4.2 The Partial Enumeration Method Lawler and Belle (1966) present a method for solving discrete optimization problems by using the concept of partial enumeration. The method is applicable to any problem with a monotonic objective function and arbitrary constraints. To be more specific, the type of problems that can be solved must have a following form:

Minimize

g0 (x)

subject to g 11 ( x ) ? g 12 ( x ) ≥ 0 g 21 ( x ) ? g 22 ( x ) ≥ 0 M g m1 ( x ) ? g m 2 ( x ) ≥ 0 x = ( x1 , x 2 ,..., x n ) where x j ∈ {0,1} , j = 1, 2, …, n

(4.1)

and each of the functions g0, gi1, and gi2 for i = 1, 2, …, m must be monotonic nondecreasing with respect to each of x1, x2,…, xn.

48

4.2.1 Partial Ordering and Numerical Ordering of Binary Vectors Since the decision variables, x1, x2,…, xn, are either 0 or 1, the vector x = (x1, x2,…, xn) is a binary vector. The lexicographic or numerical ordering of a binary vector is obtained by identifying the numerical value of each binary vector by the following formula:

n( x ) = x1 2 n ?1 + x2 2 n ? 2 + ... + xn 2 0 .

One approach for finding the optimal solution to the optimization problem is to arrange all possible solutions in numerical order and analyze them one at a time. Since the decision variables are binary numbers, there are 2n possible solutions to analyze. Lawler and Belle’s algorithm uses the concept of partial ordering to provide three rules to skip searching many possible solutions, thus reducing the number of computational steps. According to their definition of vector partial ordering, given two binary vectors x and y, y is comparable to x, denoted x ? y, if and only if xi ≤ yi, for all i = 1, 2, …, n. In addition, two binary vectors that are not comparable will be denoted as x ? y. For example, (0 1 0 0) ? (0 1 1 1) (0 0 1 1) ? (0 1 0 0) It should be noted that the numerical ordering is a refinement of the partial ordering, i.e., x ? y implies n(x) ≤ n(y): however, n(x) ≤ n(y) does not imply x ? y. For examples, n(0 1 0 0) ≤ n(0 1 1 1) and (0 1 0 0) ? (0 1 1 1) n(0 0 1 1) ≤ n(0 1 0 0) but (0 0 1 1) ? (0 1 0 0). 49

4.2.2 Skipping Rules The Lawler and Belle algorithm is an adaptation of total enumeration in numerical increasing order of binary vectors. Suppose all possible solutions are arranged in numerical order from the smallest to the largest numerical value, i.e. (0, . . ., 0, 0, 0, 0) (0, . . ., 0, 0, 0, 1) (0, . . ., 0, 0, 1, 0) (0, . . ., 0, 0, 1, 1) : (1, . . ., 1, 1, 1, 1)

It should be noted that the vectors following an arbitrary vector x in the numerical ordering may or may not be comparable to x. For any given x, treat x as a binary number when using addition (+) and subtraction (-) operations. Let x* be the first noncomparable vector following x in numerical ordering, then the vector x* can be found as follows: 1. Subtract 1 from x to obtain x-1, 2. logically ‘or’ x and x-1 to obtain x*-1, 3. Add 1 to x*-1 to obtain x*. Note that, since x* is the first vector that is not comparable to x, then x+1, x+2, …., x*-1 are all comparable to x, and x*-1 is comparable to x, x+1, x+2, …., x*-2.

50

To illustrate the skipping rules of Lawler and Belle’s method, consider a simple optimization problem: Minimize g(x) subject to ? a( x ) ≥ 0 b( x ) ≥ 0 x = ( x1 , x2 ,..., xn ) x j ∈ {0,1} for j = 1,2, ..., n Lawler and Belle’s method searches from the smallest numerical value to the largest numerical value, keeping a record of the best feasible solution found so far. Three rules listed below indicate conditions that permit skipping from the current vector x to x*. Let x be the current solution, and x+ be the best feasible solution found so far. The skipping rules are: (4.2)

1) If g(x) ≥ g(x+), skip to x*. Because g is monotonic nondecreasing, and x is not better than x+, then none of the x+1, x+2,…, x*-1 can be better than x. For example, consider the numerical ordering of these vectors.

Current x (0 0 1 0 0) g(x) x+1 (0 0 1 0 1) g(x) ≤ g(x+1) x+2 (0 0 1 1 0) g(x) ≤ g(x+2) x*-1 (0 0 11 1) g(x) ≤ g(x*-1) x* (0 1 0 0 0) ?

51

2) If g(x) < g(x+) and x is a feasible solution, skip to x*. Because g is monotonic nondecreasing, none of the x+1, x+2,…, x*-1 can be better than x. Since g(x) < g(x+) and x is feasible, x is substituted for x+. 3) If g(x) < g(x+) and –a(x)<0 or b(x*-1) < 0, skip to x*. Because a is monotonic nondecreasing, then –a is monotonic nonincreasing, if -a(x) < 0 , then none of x+1, x+2,..,x*-1 can be greater than or equal to zero. For the same reason, b is monotonic nondecreasing, if b(x*-1) < 0, then none of x+1, x+2,.., x*-1 can be greater than zero. For example,

Current x (0 0 1 0 0) -a(x) b(x) ≤ b(x*-1) x+1 (0 0 1 0 1) -a(x) ≥ -a(x+1) b(x+1) ≤ b(x*-1) x+2 (0 0 1 1 0) -a(x) ≥ -a(x+2) b(x+2) ≤ b(x*-1) x*-1 (0 0 11 1) -a(x) ≥ -a(x*-1) b(x*-1) x* (0 1 0 0 0) ? ?

If none of the three rules applied, then replace x with x+1 (no skip).

Lawler and Belle also suggest that the algorithm is more efficient if those variables that are ‘least significant’ are assigned positions at the right-hand end of the binary vector.

52

4.3 Model Preparation Before applying Lawler and Belle’s algorithm to the explosive screening device allocation model, it must be consistent with the form (4.1) or the simpler form (4.2). The model can be transformed to standard form (4.2) by the following steps:

i.

Change in decision variables The model has both integer (number of screening devices required for each type)

and binary variables. To eliminate the process of selecting the number of devices required, Yj, recall capacity constraints (3.3) and (3.4)

C j Y j ≥ ∑ Bi X ij , for j = 1, 2,…, D,

i =1 M

C j Y j < ∑ Bi X ij + C j , for j = 1, 2,…, D,

i =1

M

.

Yj ∈ Ζ + .

Let Y'j be a real-valued variable greater than or equal to zero, substitute for Yj in both equations, we know that Y'j where, Y ' j C j = ∑ X ij Bi , satisfies both constraints:

i =1 M

C j Y ' j ≥ ∑ Bi X ij , for j = 1, 2,…, D,

i =1

M

C j Y ' j < ∑ Bi X ij + C j , for j = 1, 2,…, D.

i =1

M

53

Then we have ?M ? Y ' j = ? ∑ X ij Bi ? C j , for j = 1, 2,…, D. ? i =1 ? The number of devices, Yj, must be integer, and can be easily found by rounding Y'j up to the nearest integer value. This is done by adding the following two constraints, ?M ? Yj ≥ Y ' j = ? ∑ X ij Bi ? C j , ? i =1 ? ?M ? Yj < 1+ Y'j = 1 + ? ∑ X ij Bi ? C j . ? i =1 ? The explosive screening device allocation problem now becomes maximize subject to

D

∑L R

i =1 i

M

i

∑ (Y F

j =1 j

j

+ S j I j ) + ∑∑ X ij O j B j ≤ TB

M D i =1 j =1

Where Y'j ≤ Yj < Y'j +1, j = 1,2, . . ., D Xij ∈ {0,1} Yj ∈ Ζ + . D 1 and Li = 1 ? ∏ (1 ? X ij Pj ) Ri = Ai Bi M j =1 ∑ Ai Bi

i =1

(4.3)

Y'j =

∑X

i =1

M

ij

Bi

, j = 1,2, . . ., D

Cj

This method makes Yj to be dependent variable (depends on Xij), by finding Y'j , and rounding up the real value variable, Y'j, to the integer variable Yj. Therefore, for each

54

combination of Xij, Yj can be obtained directly, thus, Yj is no longer a decision variable. Now the problem is reduced to a Binary Nonlinear Program.

ii. Change in objective function and constraints The objective function must be in the form of minimizing a nondecreasing function. Since, the model objective function is

D M ? ? Maximize ∑ ?1 ? ∏ (1 ? X ij Pj )?Ri , i =1 ? j =1 ?

it can be changed to a minimization o by multiplying by -1, yielding

D M ? ? Minimize ? ∑ ?1 ? ∏ (1 ? X ij Pj )?Ri . i =1 ? j =1 ?

(4.4)

However, the objective function is also converted into a nonincreasing function. To change it to a nondecreasing function, let X'ij = 1- Xij, and substitute into (4.4). The objective function becomes

D M ? ? Minimize ? ∑ ?1 ? ∏ (1 ? (1 ? X 'ij ) Pj )?Ri i =1 ? j =1 ?

(4.5)

which minimizes a nondecreasing function as desired. To change the budget constraint to be greater than or equal to zero, multiply both sides by -1 and move –TB to the left-hand side to obtain

TB ? ∑ (Y j F j + S j I j ) ? ∑∑ (1 ? X 'ij )O j B j ≥ 0 .

D M D j =1 i =1 j =1

(4.6)

55

Next, replace the objective function and constraint with Equation (4.5) and (4.6). The problem becomes

D M ? ? Minimize ? ∑ ?1 ? ∏ (1 ? (1 ? X 'ij ) Pj )?Ri i =1 ? j =1 ?

subject to

TB ? ∑ (Y j F j + S j I j ) ? ∑∑ (1 ? X 'ij )O j B j ≥ 0

D M D j =1 i =1 j =1

Where Y'j ≤ Yj < Y'j +1 Y'j ≥ 0 Yj ∈ Ζ + Xij ∈ {0,1} X'ij = 1- Xij 1 , i = 1, 2,…, M Ri = Ai Bi M ∑ Ai Bi

i =1

(4.7)

Y'j =

∑ (1 ? X '

i =1

M

ij

) Bi

Cj

, j = 1, 2, . . ., D

The problem is now in the standard form of (4.2) and can be solved by Lawler and Belle’s algorithm with

D M ? ? g ( x) = ?∑ ?1 ? ∏ (1 ? (1 ? X 'ij ) Pj )?Ri , i =1 ? j =1 ?

b( x) = TB ? ∑ (Y j F j + S j I j ) ? ∑∑ (1 ? X 'ij )O j B j ,

D M D j =1 i =1 j =1

a ( x) = 0 .

56

It should be noted that the equation

b( x) = TB ? ∑ (Y j F j + S j I j ) ? ∑∑ (1 ? X 'ij )O j B j

D M D j =1 i =1 j =1

is monotonic nondecreasing because each term is monotonic nondecreasing. To show this, note that the second term,

M D

∑∑ X

i =1 j =1

M

D

ij

O j B j is monotonic nondecreasing in each Xij,

then

M

∑∑ (1 ? X '

i =1 j =1 D

ij

)O j B j is monotonic nonincreasing in each X'ij, as a result,

? ∑∑ (1 ? X ' ij )O j B j is monotonic nondecreasing in each X'ij.

i =1 j =1

Also for the term ? ∑ (Y j F j + S j I j ) , consider the equation

D j =1

Y'j =

∑X

i =1

M

ij

Bi

.

Cj

Note that Y'j is nondecreasing if Xij is increasing. Then

Y'j =

∑ (1 ? X '

i =1

M

ij

) Bi

Cj

is monotonic nonincreasing for each X'ij. Consequently,

D

∑ (Y F

D j =1 j

j

+ S j I j ) is monotonic

nonincreasing for each X'ij. As a result ? ∑ (Y j F j + S j I j ) is monotonic nondecreasing for

j =1

each X'ij.

57

4.4 Data Generation In this research, the airport data and the screening device data were determined using information available in the existing literature (Butler and Poole, 2002; McLay, 2006; Virta et al., 2003) and experimental justification. Some information may not be available as it is classified. Therefore, some parameters which are not available are computed based on the available information and estimation. For the screening device, the fixed and installation costs are obtained based on the expected useful life of the device and the number of operational hours per year (per day times 365). These costs are computed by the yearly costs divided by hours of operation. In brief, the computational steps for randomly generating the screening device data are as follows. 1. Probability of threat detection is randomly assigned from a uniform distribution ranges from 0.5 to 0.99. 2. Capacity is normally distributed with mean and standard deviation equal to 100 and 50, respectively. 3. Fixed cost is 20% of capacity times the uniform random number ranges from 0.5 to 1.5. 4. Operating cost is approximately 5% of fixed cost. 5. Installation cost is approximately 4 times the fixed cost.

The underlying ideas for generating airport data are roughly outlined as follows. 1. Airport budget is randomly generated from a triangular distribution with minimum, mean, and max equal to 500, 5000, 1200, respectively.

58

2. Total number of passengers is computed by dividing the budget by a uniformly distributed integer random number from 2 to 4. 3. The number of classes is assigned with probability 0.9 for having 2 classes, and 0.1 for having 3 classes. 4. The number of passengers in each class is computed as, ? For a 2-class airport, approximately 85%, and 15% of total passengers are assigned to class 1 and class 2, respectively. ? For a 3-class airport, approximately 70%, 20%, and 10% of total passengers are assigned to class 1, 2, and 3, respectively. 5. The average assessed threat value for class 1, 2, and 3, are independent and normally distributed with means equal to 0.23, 0.35, and 0.55, respectively. 6. For the existing number of devices, the probability that an airport has each device is equal to 0.3, then the number of devices is computed as the total number of passengers divided by the device capacity.

There is an additional requirement for passenger class: The average assessed threat value of class 1 must be strictly less than class 2, and class 2 must also be strictly less than class 3. All numerical parameters listed above may vary based on experimental justification and appropriateness of data. The randomly generated data for both airports and screening devices are used throughout the following section and in Chapter 7.

59

4.5 Computational Examples In this section, computational testing is used to evaluate the performance of the algorithm. The method is demonstrated by various sizes of examples (number of devices and number of classes). The task is to find a combination of explosive screening devices for each class in order to maximize the total security. Representative data for seven screening devices and four classes are randomly generated and given in Table 4.1 and 4.2, respectively.

Table 4.1. Screening Devices Data. Probability Device of threat detection 1 2 3 4 5 6 7 0.80 0.70 0.60 0.75 0.95 0.68 0.78 Number Capacity (bag/device) 150 28 200 40 100 200 120 Operating cost ($/bag) 1.00 0.83 0.69 0.84 0.94 0.50 0.70 Fixed Install of devices 41 11 10 20 60 2 30 100 40 15 70 250 20 30 1 0 2 1 2 0 1 ($/device) ($/device) existing

60

Table 4.2. Passenger Classes Data. Average assessed threat value 0.20 0.44 0.70 0.10

Class 1 2 3 4

Number of bags 350 150 80 500

The Lawler and Belle’s algorithm was implemented in MATLAB, and was tested on an Intel Pentium 4 processor, 2.4 GHz with 512 MB of RAM. The method starts by assuming the best solution found so far, x+, is (0 0 …. 0 0), the best objective value is g(x+) = - K, where K is a very large number, and the first binary vector to consider is (0 0 0 … 0), where the first D (number of device types) components of the vector represent X11 to X1D, the second D components represent X21 to X2D, and so on. The three rules are tested against the current solution until the algorithm is stopped when the current solution is 1(0 0 0 … 0), namely overflow (or the vector (1 1 1 … 1) + 1). A flow chart of the computational steps is shown in Fig 4.1. Computational results for various sizes of the problems are shown in Table 4.3. The first two columns show the problem size, e.g. for 2 x 2, the data is from the first two classes of passengers and the first two devices from Table 4.2 and 4.1, respectively. For convenience, solutions are shown in matrix form [Xij]MxD, i=1, 2, …, M, and j = 1, 2, . . ., D, where Xij represents the ith class of passengers being screened by device type j.

61

START

x = (0 0 …0) x+=1(0 0 …0) g(x+)= - K

Set x = x+1

x = 1(0 0 …0)? (over flow) Yes No STOP

Set x = x*

g(x) ≥ g(x+)? (compare obj)

Yes (rule 1)

No

b(x) ≥ 0? (feasibility)

Yes (rule 2)

Set x+= x

No

No

b(x*-1) ≥ 0? (feasibility)

Yes (rule 3)

Figure 4.1. Flowchart of Lawler and Belle’s Algorithm.

62

Table 4.3. Computational Results. Problem size # of classes (M) # of devices (D) 2 2 Budget 2000 Solution 1 1 1 0 Y = [4 13] S = [3 13] 1 0 1 1 1 1 Y = [4 6 3] S = [3 6 1] 1 0 1 1 1 1 1 0 1 Y = [4 6 3] S = [3 6 1] 0 0 0 0 1 1 0 0 0 1 1 0 1 1 0 Y = [2 0 1 2 5] S = [1 0 0 1 3] 1 0 0 0 0 1 1 0 0 0 0 1 1 1 1 0 1 1 0 1 1 Y = [3 0 1 2 2 3 5] S = [2 0 0 1 0 3 4] 1 0 1 0 1 0 1 0 1 0 1 1 1 0 1 0 Y = [8 0 6 2] S = [7 0 4 1] 1 0 1 0 0 1 0 0 1 0 1 1 1 0 1 0 0 1 1 0 0 0 0 1 Y = [7 0 3 0 2 6] S = [6 0 1 0 0 6] Objective value 0.87206 Total cost 1933.05

2

3

2000

0.94718

1793.12

3

3

2000

0.93925

1928.32

3

5

2200

0.9725

2175.79

3

7

2300

0.99325

2290.41

4

4

3200

0.93388

3155.76

4

6

3200

0.9717

3199.89

63

The results of Lawler and Belle’s algorithm were verified by the Search All method, which enumerates all combinations from the lowest numerical value, i.e. (0 0 0 … 0), to the highest numerical value, (1 1 1 . . .1) and keeps the best solution. The comparison of the two algorithms, which includes size of the problem, number of possible solutions, number of iterations using the LB algorithm, and computational time, are given in Table 4.4.

Table 4.4. Comparison of Two Algorithms. Total Problem size Number of classes (M) 2 2 2 3 3 3 3 4 4 Number of devices (D) 2 3 4 3 5 6 7 5 6 16 64 256 512 32,768 262,144 2,097,152 1,048,576 16,777,216 =2(MxD) possible solutions Lawler and Bell’s algorithm # of Iterations 6 12 74 44 1,591 3,463 6,546 70,328 489,664 LB 0.031 0.047 0.062 0.047 1.329 2.938 5.719 58.469 429.21 Search All 0.016 0.047 0.172 0.266 15.593 130.562 1003.210 503.872 7980.425 YES YES YES YES YES YES YES YES YES Elapsed time(sec) Optimal ?

64

The first two columns of Table 4.4 show the problem size in number of classes and number of devices. The column labeled Total Possible Solutions is the number of possible combinations of (M x D) binary numbers, which is equal to the number of Search All steps. The next column shows the number of iterations of Lawler and Belle’s algorithm, i.e. the number of LB steps. Elapsed Time is the running time of the computer for each algorithm in seconds. The last column checks whether both algorithm solutions are the same. It can be seen that the running times of both algorithms increase exponentially as the size of the problem increases. Currently, passengers are classified into only two classes, and there are less than ten types of baggage screening devices to choose from. For the problem size of 4 x 6, the Lawler and Belle method takes less than 500 seconds (less than 10 minutes) to obtain the optimal solution. Therefore, the computational results show that the Lawler and Belle method is a robust and easy to implement method for solving a reasonably sized singleairport problem.

65

CHAPTER VI MULTIPLE AIRPORT METHODOLOGY

6.1 Introduction In this chapter, problem solving approaches for multiple airports model are presented. Grouping many airports and attempting to simultaneously solve the set of all airports complicate the problem because the number of binary decision variables is dramatically increased. The model is NP-hard and, from a practical standpoint, is extremely difficult to solve to optimality, even for a relatively small size problem. For example, suppose there are 30 airports under consideration, 4 types of devices, and 2 classes of passengers, then there will be 30 x 4 x 2 = 240 binary decision variables (not including the number of screening device decision variables), therefore there are approximately 2240 ≈ 1.66 x 1072 possible configurations of the devices. The optimization problem becomes a large-scale binary problem with a nonlinear objective function and multiple constraints. The Lawler and Belle algorithm obviously cannot handle this problem within reasonable time. New algorithms or heuristics which are less timeconsuming must be studied to solve this large scale problem. The majority of the screening device and passenger data are estimated using the numbers that represent the best educated guess about their expected value. They usually contain errors of approximation. The probability of threat detection and the passenger’s assessed threat value are both subject to random fluctuation. An optimal solution for the estimated data may not actually be an optimal solution for the real data. Solving the

74

model by an exact algorithm to optimality with approximate data may not be worthwhile in this situation. Therefore, we focus on developing heuristics instead of an exact algorithm. To find an acceptable solution for the multiple airports problem, three heuristic approaches are developed in this chapter. The first heuristic approach is based on a greedy heuristic in which the airport with the highest number of bags is iteratively chosen to be solved to optimality. The second approach focuses on greedily choosing each binary decision variable by considering the sorting and updating of sensitivity factors. The third approach is based on creating and solving a multidimensional knapsack problem combined with the first or second approach.

6.2 The Effect of the Interaction Variables It is obvious that the overall objective function of the Multiple Airport Model,

∑ SLk , where SLk = ∑ Bik Lik Ri k ,

k =1

T

Mk

i =1

is separable, i.e. it is the sum of each airport security level. Unfortunately, there are interaction variables in the model constraints (the number of new devices needed); optimizing each airport independently without regard to the effects of interactions may not achieve an overall optimal solution. However, there are some circumstances where the interactions may have little or no effect in the model. These cases include 1) There are no new devices available. 2) There are more than enough new devices available for all airports.

75

3) The cost of installing a new device is very high and not justified for most airports. To better explain the situations, consider the new device constraints, Equation (5.9), which are binding (connecting) constraints,

∑S

k =1

T

jk

≤ U j , for j = 1,2, . . ., D.

In case 1), there will be no interactions among airports because each Uj will be zero (no new device available) forcing each Sjk to be zero. For case 2), there are interactions but they are insignificant, since they will never violate new device constraints (the sum of Sjk will not be greater than Uj ). For case 3), the optimal solution may contain most of Sjk as zero (installing new device may be too expensive for most airports). Without the interaction variables, Sjk, the optimal solution of the overall system is achieved by optimizing each subsystem or airport. This procedure may yield a good overall solution, if the effects of interactions are insignificant (as in case 3).

6.3 Airport Greedy Approach The objective of this section is to introduce the Airport Greedy Heuristic which selects airports that will be optimized one at a time. Airport selection is based on the number of bags or passengers in the airport. We solve each airport to optimality and update the availability of the new devices, Uj, after the optimal solution for each airport is

76

found. The process is repeated until there are no airports left to consider. This heuristic will be referred as the AGH for the remaining of this dissertation. The AGH steps are described as follows. Let A denote the set of airports under consideration, Z be the overall objective function value, and Zk be the security level of airport k. Initially A= {1, 2, 3,…., T} and Z =0. ? ? Step 1: Select airport, k, in set A which has the most bags (or passengers). Step 2: Solve airport k using the Lawler and Belle algorithm with airport k parameters (budget, number of passengers, number of classes, and number of existing devices) and number of devices available Uj, to obtain Zk and Sjk, for j = 1,2, . . ., D. ? ? Step 3: Update the resource availability, Uj = Uj - Sjk, for j = 1,2, . . ., D. Step 4: Update the overall objective function value, Z = Z + Zk, and remove airport k from set A. ? Step 5: If set A is empty, stop, otherwise go to step 1.

AGH can also be described more precisely by the control abstraction below. Algorithm AGH A = {1, 2, 3, …., T} //set of airports under consideration Z=0 //current overall objective function value while A ≠ {} select k from A, where BAGS(k) = max{BAGS(i)| i ? A} [Sjk, Zk ] = OPTIMIZE({airport k parameters}, {Uj}) for j = 1 to D Uj = Uj - Sjk repeat Z = Z + Zk A = A – {k} end while end AGH Figure 6.1. Algorithm AGH Control Abstraction.

77

The function BAGS calculates the total number of bags for airport k, OPTIMIZE finds the optimal solution with airport k parameters and the current number of devices available, and returns the number of devices to install and the security level of the airport.

6.4 Variable Greedy Approach The only decision variables in the Multiple Airport model are Xijk, binary decision variables representing the decision of using device type j to screen class i bags at airport k, where Yjk, the number of device type j required at airport k, can be assigned optimally as discussed in Chapter 4. In this section, we present a greedy heuristic approach which iteratively selects a binary decision variable Xijk that has the maximum sensitivity factor. In the beginning, all decision variables are set to be zero. Then, the sensitivity factor for each decision variable is calculated. The sensitivity factor used in this algorithm is the increment in the overall objective function value per unit of capacity used. The solution is obtained by repeatedly selecting the decision variable Xijk which has the maximum sensitivity factor value without violating any of the constraints. First, we define the function OBJ(Xijk) to be the increment in objective value if Xijk is changed from zero to one. To determine the value of this function, consider Equation (5.4), the objective function for each airport,

SLk = ∑ Bik Lik Ri k ,

i =1

Mk

where Bik, Rik, and Lik, are the number of bags, the risk level, and the probability of threat detection for class i bags at airport k, respectively. As the value of Xijk is changed to one,

78

it can be seen from Chapter 3 that Bik and Rik are fixed, while Lik increases. To analyze the change in Lik, consider the Lik equation, Lik = 1 ? ∏ (1 ? X ijk Pj ).

D j =1

Initially all Xijk are set to zero, therefore Lik will also be equal to zero. It is obvious that if Xijk is changed to one, Lik is increased to Pj, therefore the increment in Lik is equal to Pj. However, if there is a binary variable Xijk that is already selected (set to one) in the same class or same airport, the increment in Lik is not so straightforwardly obtained. The value of the increase in Lik can be demonstrated by the following example. Suppose that the device type 1 of class i at the airport k, is already selected, then Lik is equal to P1. If device type 2 is additionally selected for the same class and airport, then Lik is equal to Lik after device 2 is selected = 1- (1- P1)(1- P2). Therefore the increment in Lik is equal to the probability of threat detection after selection of device 2 minus the probability before selection of device 2 or, Lik(after) - Lik(before) = [1- (1- P1)(1- P2)] - P1. Denote this increment as ?Lik, then the rate of improvement in the objective value if Xijk is selected is the product of ?Lik, Rik, and Bik, or OBJ(Xijk) = ?Lik·Rik·Bik. (6.1)

There are two types of cost in the Multiple Airport Model; the cost of using the device and installing a new device. In order to compute both types of costs, recall the cost equation

∑ (Y

D j =1

jk F j + S jk I j ) + ∑∑ X ijk O j Bik . Mk D i =1 j =1

79

The value of Yjk, the total number of devices needed, and Sjk, the number of new devices to be installed, are still undetermined. The number of devices needed can be approximated by the number of bags screened by that device type divided by the capacity of the device, denote Y`jk as the approximate number of devices needed, then Y`jk = Bik / Cj. Determine the approximate number of new devices to be installed S`jk = Y`jk – Ejk if Y`jk > Ejk, and Sjk = 0 otherwise, then the approximate cost increase of setting Xijk = 1 can be directly computed as

Y `jk F j + S `jk I j + O j Bik .

Each airport has a different budget, therefore the cost increment for each variable should be scaled in order to have the same magnitude. Divide the cost increment by the total budget of each airport, in order to obtain the cost ratio, or Cost(Xijk), Cost(Xijk) = (Y `jk F j + S `jk I j + O j Bik )/ TBk . Similarly, the new device resource consumption ratio, denoted New(Xijk), can be defined as the number of new devices to be installed divided by the total number of new devices of that type available, or New(Xijk) = S `jk / U j . It is possible that there may be no resources available for the new device. In such a case, to avoid dividing by zero, we assign New(Xijk) = 1 if S`jk > 0, and New(Xijk) = 0 if S`jk = 0. It should be noted that after the decision variable Xijk is selected, the number of existing devices (if any) of type j at airport k changes. The cost ratio and the new device ratio of 80

decision variables for the same airport, k, and device type, j, must be updated after each selection. The outline of the Variable Greedy Heuristic, referred as VGH, is as follows. Initially, set all decision variables, Xijk, and the current optimal solution value, Z, equal to zero. ? Step 1: Calculate the increment in objective function, OBJ(Xijk) = Pj·Rik·Bik, for all decision variables. ? Step 2: Calculate the approximate number of devices used, Y`jk = Bik / Cj, for all decision variables. ? Step 3: Calculate the approximate number of new devices to be installed, S`jk = Y`jk – Ejk if Y`jk > Ejk, and Sjk = 0 otherwise, for all decision variables. ? Step 4: Calculate the increment ratio of cost, Cost(Xijk) = (Y `jk F j + S `jk I j + O j Bik )/ TBk , for all decision variables. ? Step 5: Calculate the increment ratio of the new device, if Uj > 0, New(Xijk) = S ` jk / U j , if Uj = 0 and S`jk > 0, New(Xijk) = 1, if Uj = 0 and S`jk = 0 New(Xijk) = 0, for all decision variables. ? Step 6: Calculate the sensitivity factor, FAC(Xijk), as FAC(Xijk) = OBJ(Xijk) / (Cost(Xijk) + New(Xijk)), for all decision variables. ? Step 7: Select the decision variable Xijk which has the maximum sensitivity factor. If the maximum sensitivity factor is equal to zero, take current solution as the VGH solution and stop. Otherwise, continue on step 8. 81

?

Step 8: Test whether the decision variable selected in step 7 is feasible (does not violate the constraints). If it is feasible, set Xijk = 1, FAC(Xijk) = 0, Z = Z + OBJ(Xijk) and go to step 9. Otherwise, set FAC(Xijk) = 0 and go to step 7.

?

Step 9: For the same airport and class as the decision variable selected in step 7, recalculate sensitivity factors by updating the increment in the objective function, OBJ(Xijk) = ?Lik·Rik·Bik, where ?Lik = Lik(after select) - Lik(before select).

?

Step 10: For the same airport and device type as the decision variable selected in step 7, recalculate the sensitivity factors by o Updating the remaining number of existing device, Ejk = Ejk – Y`jk if Ejk > Y`jk, otherwise Ejk =0. o Then, update the cost ratio, new device ratio and sensitivity factor with the new Ejk, and go to step 7. The VGH algorithm is best described by the control abstraction in Figure 6.2.

Function Update_Objective updates the increase in the objective values, OBJ(Xijk), for all device types at the same airport and class of previously selected variable. Function Update_Ex_Device determines the number remaining devices, Ejk, for all classes at the same airport and device type of the previously selected variable. Function Update_Cost and Update_Cost_Ratio update the cost increment and cost ratio, Cost(Xijk), respectively, of all classes at the same airport and device type of previously selected variable. Function Update_Factor updates the factor with the new OBJ(Xijk) and Cost(Xijk) for all classes and device types at the same airport of the previously selected variable.

82

Algorithm VGH A = {1, 2, 3, …., T} //set of airports under consideration Z =0 //current overall objective function value Mk = the number of passenger classes at airport k, k ? A. D = the number of screening device types. X = { Xijk } // set of current solution for k =1:T, j=1:D, i =1:Mk OBJ(Xijk) = Pj·Rik·Bik //increment in objective Y`jk = Bik / Cj //approx. number of device used if Y`jk > Ejk, S`jk = Y`jk – Ejk else Sjk = 0 end if Cost(Xijk) = (Y `jk F j + S `jk I j + O j Bik )/ TBk //cost increment ratio if Uj > 0 New(Xijk) = S `jk / U j else if Uj = 0 and S`jk > 0 New(Xijk) = 1 else if Uj = 0 and S`jk = 0 New(Xijk) = 0 end if FAC(Xijk) = OBJ(Xijk) / (Cost(Xijk) + New(Xijk)) repeat while max{FAC (Xijk)| k =1:T, j=1:D, i =1:Mk } ≠ 0 select Xabc, where FAC(Xabc) = max{FAC(Xijk)| k =1:T, j=1:D, i =1:Mk} if FEASIBLE(Xabc) Xabc = 1 Z = Z + OBJ(Xijk) Update_Objective(class = a, device = ALL, airport = c) Update_Ex_Device(class = ALL, device = b, airport = c) Update_Cost (class = ALL, device = b, airport = c) Update_Cost_Ratio(class = ALL, device = b, airport = c) Update_Factor(class = ALL, device = ALL, airport = c) else Xabc = 0 end if FAC(Xabc) = 0 end while end VGH

Figure 6.2. Algorithm VGH Control Abstraction.

83

6.5 3-Phase Approach In Section 6.2, we heuristically pick airports then optimally solve them one at a time while updating the resource availability after each iteration. In this section, instead of selecting one airport at a time, we introduce an approximate algorithm which chooses a set of airports to be optimized without violating the new screening device resource constraints. The algorithm heuristically chooses a set of airports from the set of airports under consideration. The underlying idea of this algorithm is that the set of airports to be optimized are optimally chosen such that the solution does not violate the budget constraint for each airport or the new device resource constraints. The term “optimally chosen” means that there are no other sets that give a better objective function value and are feasible. The algorithm, referred to as the 3-Phase algorithm, is inspired by the previous work of Chern and Jan (1986). It can be divided into three major phases and is explained as follows. Phase I In Phase I, we decompose the T-airport problem into 2 types of T sub-problems (T is the total number of airports under consideration), denote problem type I as BPk (base problem), and type II as IPk (improve problem) where k = 1, 2, 3…, T. The base and improve problems have the same characteristics as the single airport model discussed in Chapter 3, where the difference is the new device resource constraints. For each base problem, it is assumed that there are no new devices available to install, while each improve problem will have the number of new devices available equal to total number of new devices available. Both problems are formulated as follows.

84

Problem BPk max SLk = ∑ Bik Lik Ri k

i =1 Mk

subject to

∑ (Y

D j =1

jk F j + S jk I j ) + ∑∑ X ijk O j Bik ≤ TBk , Mk D i =1 j =1 G

Y jk C j ≥ ∑ (X ijk Bik ), for j = 1,2,..., D,

i =1 M

Y jk C j < ∑ (X ijk Bik ) + C j , for j = 1,2,..., D,

i =1

S jk = 0, for j = 1,2,..., D. (or Y jk ≤ E jk , for j = 1,2,..., D) Problem IPk max SLk = ∑ Bik Lik Ri k

i =1 Mk

subject to

∑ (Y

D j =1

jk

F j + S jk I j ) + ∑∑ X ijk O j Bik ≤ TBk ,

Mk D i =1 j =1 G

Y jk C j ≥ ∑ (X ijk Bik ), for j = 1,2,..., D,

i =1 M

Y jk C j < ∑ (X ijk Bik ) + C j , for j = 1,2,..., D,

i =1

S jk ≤ U j , for j = 1,2,..., D. where

X ijk ∈ {0,1}, Y jk ∈ Z + , Lik = 1 ? ∏ (1 ? X ijk Pj ), and Rik = Aik Bik

D j =1

1

∑A

i =1

Mk

. Bik

ik

85

Both types of problems are solved to optimality with their own budget constraints, so that the optimal solution of each airport does not violate any constraints. The difference in objective values of both types of problems for each airport is calculated and referred as Profit(k). Phase II In this phase, we construct a knapsack problem with a collection of different new device resource constraints, or a multidimensional knapsack problem (D-KP), where D is the number of device types. The problem D-KP can be defined as follows.

max

∑ Profit(k ) ? u

k =1 jk

T

k

subject to

∑S

k =1

T

? u k ≤ U j , j = 1,..., D,

u k ∈ {0,1}, k = 1,..., T

where Sjk is the new device type j installed at airport k from the improve problem, Profit(k) is the difference in the objective value of IPk and BPk, Uj is the number of device type j available, T is the total number of airports under consideration, and D is the total number of new device types. The binary decision variable, uk, represents the decision about using the optimal solution from the improve problem. We solve the 0-1 multidimensional knapsack problem and obtain the optimal solution vector u = { u1, u2, …, uk}. Phase III In Phase III, for the set of airports that have solution uk = 1, denoted the set of optimized airports, the airport solutions are the same as the solution obtained from solving the improve problem in Phase I. There might be the case where there are plenty

86

of new devices available (as case 3 mentioned in Section 6.2). In this case, the sum of all the weights is less than or equal to the resource available, or

∑S

k =1

T

jk

≤ U j , for all j, then

clearly uk = 1, for all k, is an optimal solution, and the algorithm is terminated early. For the set of airports with uk = 0, denoted the set of unoptimized airports, two cases are considered. Case one, if there are no new device resources left, the airport solutions are obtained from the base airport in Phase I. However, in Case two, if there are new devices remaining; we solve the set of unoptimized airports with the remaining new device resources using the algorithm presented in Section 6.2, referred as 3PH-A, or using the algorithm in Section 6.3, referred as 3PH-V. The steps of the 3-Phase algorithm are outlined as follows. ? Step 1: Solve each improve problem to optimality with the new device resource constraints, assign the optimal objective value to Z(IPk). ? Step 2: From new resource consumption Sjk obtained from Step 1, test whether the current solution is optimal. If ∑ S jk ≤ U j , then the current solution is optimal and

k =1 T

stop, otherwise continue to Step 3. ? Step 3: Solve each base problem to optimality with the assumption that there are no new devices available to install, i.e., Uj = 0. Assign the optimal objective value to Z(BPk). ? Step 4: Calculate the difference in an objective function between IPk and BPk as Profit(k) = Z(IPk) - Z(BPk).

87

?

Step 5: Generate the multidimensional knapsack problem (D-KP), with T binary decision variables (items or airports), and D constraints (device type), where profit for item k is Profit(k), and resource consumption of item k for constraint j is Sjk (from improve airport solution).

? ?

Step 6: Solve D-KP generated from Step 5, obtain the solution vector u. Step 7: Assign solutions from improve airport problems to the airport that has binary decision variable =1. If U j ? ∑ S jk ? u k = 0 , for all j, assign solutions from

k =1 T

base airport problems for the rest of the airports (uk = 0), and stop. Otherwise continue to Step 8. ? Step 8: Solve the set of airports that have uk = 0, by algorithm VGH (or AGH) with updated resource constraints, U j ? ∑ S jk ? u k , then stop.

k =1 T

Algorithm 3PH can be explained by the control abstraction in Figure 6.3. Method OPTIMIZE_1 finds the optimal solution of the single airport problem, keeps the solution, and returns the objective value and number of new devices allocated. OPTIMIZE_2 performs the same functions as OPTIMIZE_1, except that it does not return the number of new devices allocated. Method SOLVE_KNAPSACK finds an optimal solution for the multidimensional knapsack problem and returns the optimal solution, where the profit of each item is the profit of choosing the airport optimal solution, the weight of each item is the number new devices allocated obtained from method OPTIMIZE_1, and the dimension of the constraints is the number of device types. The function UNION adds an airport that has uk = 0 to an unoptimized airports set, UN. Variable resource_remain(j) is 88

the number of new devices remaining after assigning new devices to a set of optimized airports.

Algorithm 3PH A = {1, 2, 3, …., T} // set of airports under consideration Z=0 u = {0} // D-KP solution vector xk = {0} // airport k solution vector UN = {} // set of unoptimized airports for k = 1:T [Sjk, Z(IPk) ] = OPTIMIZE_1({airport k parameters}, {Uj}) repeat if

∑S

k =1

T

jk

≤U j

u = {1, 1, 1, …,1} STOP 3PH end if for k = 1:T Z(BPk) = OPTIMIZE_2({airport k parameter}, {Uj=0}) Profit(k) = Z(IPk) - Z(BPk) repeat u = SOLVE_KNAPSACK(profit={Profit}, weight={Sjk},B={Uj}) for k = 1:T if uk = 1 xk = SOLUTION(Z(IPk)) else UN = UNION(UN, uk) end if repeat for j =1:D

resource_remain(j) = U j ? ∑ S jk u k

k =1

T

repeat call VGH (A = UN, Uj = resource remain(j), for j=1:D) (call AGH(A = UN, Uj = resource remain(j), for j=1:D)) end 3PH

Figure 6.3. Algorithm 3PH Control Abstraction. 89

6.6 Chapter Summary In this chapter, problem solving schema as well as three heuristics are presented. The first heuristic, AGH, is based on greedily selecting each airport to optimize. The essence of the second heuristic, VGH, is greedily selecting each variable using the assigned sensitivity factor. For the last heuristic, 3PH, we construct a multidimensional knapsack problem (D-KP), and solve it to optimality. The binary variable in D-KP represents the decision of selecting an airport to be optimized first. The rest of the airports are grouped together and solved by the first or the second heuristic. To examine the algorithm performance in terms of both solution quality and computational time, an experiment must be conducted. The results of the experiment are discussed in the next chapter.

90

CHAPTER VII EXPERIMENTAL RESULTS FOR MULTIPLE AIRPORT PROBLEM

7.1 Introduction In Chapter 6, we presented three different approaches (four heuristics) for finding an approximate solution of the multiple-airport problem. The algorithms were coded with MATLAB (Appendix A) in order to evaluate the performance of the proposed algorithms. In this chapter, the performance of the algorithms is compared using the objective value and computation time for different test problems. For large problems, optimal solutions cannot be obtained within reasonable time by using Lawler and Belle’s algorithm. For these problems, the solution quality is measured by the percentage deviation of the objective value from the upper bound instead of the optimal value.

7.2 Exact Algorithm In order to obtain the objective value for the optimal solution, two commercial programs, What’s Best 8.0 and Premium Solver Platform for Excel 6.5, which claimed to be able to solve mixed integer nonlinear problem were tested with a multiple-airport model. Both programs are added-in to the standard Solver in Microsoft Excel. The test problem consists of two airports, each of which has two classes of passengers and three types of device (total of eight binary decision variables). The problem was created in

91

Microsoft Excel spreadsheet. Both programs found an integer feasible solution within a minute; however, they were unable to find the global optimal solution within one hour. The one airport problem is NP-Hard and normally requires a lot of computational time to solve. As the number of airports increases, the model becomes a large scale nonlinear integer programming problem, and the computational time increases exponentially. Contrary to linear programming, which is convex, integer programming problems have many local optimal solutions. To find the global optimal solution for integer programming, we must prove that the particular solution dominates all other feasible solutions. In summary, the major difficulties of the model are listed as follows. ? ? ? ? It is a large scale non-separable optimization problem. Variables are restricted to be integer (binary). The objective function is nonlinear. It contains discontinuous (piecewise linear) constraints. From the survey of existing literature, a large number of algorithms for nonlinear integer programming problem are heuristics, while the exact algorithms are computationally complicated and usually require particular characteristics of the problem, such as, convexity, twice differentiability, linear objective function, separability, etc. As discussed in Chapter 4, Lawler and Belle’s algorithm is able to solve a problem regardless of those particular characteristics. Therefore, in this chapter, the algorithm presented in Chapter 4 is revised for finding the optimal solution of the multiple-airport problem. The objective function is changed from the single-airport model to the multipleairport model. A feasibility test for new resource constraints is added to the algorithm,

92

and the budget constraint is modified to cover more than one airport. The algorithm was coded in MATLAB and the optimal solution was verified by the total enumeration method. Test problems were executed on a personal computer with an Intel Pentium 4 processor, running at a clock speed of 3.0 GHz and 2.50 GB of RAM. The computational times are the CPU time to find to the optimal solution and are shown in Table 7.1.

Table 7.1. Computational Time of Lawler and Belle’s Algorithm. Number Number Number of airports 2 3 4 5 6 of classes 2 2 2 2 2 of devices 3 3 3 3 3 Total binary variables 12 18 24 30 36 Total possible solutions 4096 262,144 16,777,216 1,073,741,824 68,719,476,736 Computational time 2 seconds 58 seconds 1334 seconds 6120 seconds 34 hours

It is obvious that, for explosive screening device allocation for the multipleairport problem, Lawler and Belle’s algorithm is faster than What’s Best 8.0 and Premium Solver Platform for Excel 6.5. Therefore, the Lawler and Belle method is used to find the optimal objective values in this dissertation.

7.3 Upper Bound It can be seen from Table 7.1 that the number of possible solutions increases exponentially in terms of the number of binary decision variables. The computational

93

time also increases dramatically as the number of airports increases. For large problems, optimal solutions for the multiple-airport problem cannot be obtained using Lawler and Belle’s algorithm within a reasonable time. The performance of the proposed methodologies is then determined by comparing its objective value to the upper bound of the optimal solution. The upper bound is derived by constructing K single-airport subproblems from a K-airport problem, then finding the upper bound of each single-airport sub-problem. Consequently, the upper bound of the multiple-airport problem is the sum of the separate single-airport optimal solution. The single-airport problem is formulated in the same manner as the Improve Problem (IPk) discussed in Section 6.4. The test problem instance is randomly generated and is shown in Table 7.2 and Table 7.3, where column Bi and ATi are the number of bags and average assessed threat value in group i, respectively, and column Ej is the number of existing devices of type j.

Table 7.2. Three Screening Device Types Data. Device 1 2 3 Probability Capacity of threat (bag/device) detection 0.80 150 0.70 28 0.60 200 Operating cost ($/bag) 1.00 0.83 0.69 Fixed Install ($/device) ($/device) 41 11 10 100 40 15

Table 7.3. Airports Data. Airport 1 2 3 4 5 Name AAA AAB AAC AAD AAE Budget 5000 4000 4600 7800 1500 B1 400 800 2074 2527 416 B2 200 300 226 73 84 AT1 0.44 0.8 0.39 0.14 0.33 AT2 0.2 0.37 0.56 0.4 0.7 E1 2 2 6 0 0 E2 4 3 26 31 0 E3 2 2 7 0 1

94

From Section 6.2, we know that there are two cases where the upper bound will be equal to the optimal solution because each airport sub-problem will attain its optimal values. These two cases are 1) if there are more than enough new devices available for all airports, and 2) if there are no new devices available at all. These two cases will be referred as two extreme cases for the rest of the chapter. Therefore, we will consider test problems with various sets of numbers of new devices available. First, we find the maximum number of new devices needed for all airports (referred as total demand) by assuming that the number of new devices available is unlimited, solving each airport independently, then finding the total number of new devices needed. The numbers of new devices available for five problems is 0%, 25%, 50%, 75%, and 100% of the total demand of new devices. The upper bound and the optimal solution for a 5-airport problem, along with the computational time are shown in Table 7.4 below.

Table 7.4. Comparison of Upper Bounds and Optimal Solutions for 5 Airports. Percentage Problem of the total demand 1 2 3 4 5 0% 25% 50% 75% 100% # of new devices available [0 0 0 ] [10 6 5] [20 13 10] [30 20 15] [40 26 20] Upper bound Optimal solution

Time (s) Objective value Time (s) Objective value 1.71 3.08 2.12 2.16 1.39 314 2765 4651 4933 4933 80 1221 2816 1085 6120 314 2440 4025 4781 4933

95

With an assumption that the number of new devices available is unlimited, the total demand is computed to be [40 26 20], i.e., 40 type 1 devices, 26 type 2 devices, and 20 type 3 devices are needed. Five test problems with different percentages of the total demand ranging from 0%-100% are used in the experiment as shown in Table 7.4. It can be seen that the computational times for the upper bounds are relatively small compared to the computational time of the exact algorithm. The quality of the upper bound is evaluated by the percentage of the relative error, calculated as Percentage Error =

UB ? OP × 100% , OP

where UB is the upper bound of the objective value and OP is the optimal solution objective value from the exact algorithm. The percentage error is shown in Table 7.5 and is plotted in Figure 7.1.

Table 7.5. Percentage Error of Upper Bound. Percentage of the total demand 0% 25% 50% 75% 100% Percentage error of upper bound 0.00% 13.32% 15.55% 3.18% 0.00%

96

40% 35%

Percentage Error

30% 25% 20% 15% 10% 5% 0% 0% 25% 50% Percentage of total demand 75% 100%

Figure 7.1. Percentage Error of Upper Bound.

It can be seen that the upper bounds obtained from the test problems are all within 20% of optimality. For the 0% and 100% problems (two extreme cases), the upper bounds are equal to the optimal values.

7.4 Experimental Result In this section, we evaluate the performance of the four heuristics presented in Chapter 6. For small size problems, 5 airports, the heuristics are compared with the optimal solution. For large problems, 100 airports, they are compared using the upper bound for the optimal solution.

97

7.4.1 Small Size Problem The test problem instances used in this section are the same as Section 7.3, the result of the experiment is shown in Table 7.6, where “Time” is the computational time in seconds, and “Obj val” is an abbreviation for the objective value obtained by each method. Table 7.7 and Figure 7.2 depict the percentage error of the four heuristics.

Table 7.6. Computational Time and Objective Value of Four Heuristics for Five-Airport Problem. % of the total demand 0% 25% 50% 75% 100% (s) 0.32 0.37 0.33 0.31 0.26 AGH Time Obj val 314 2257 3780 4431 4933 VGH Time (s) 0.41 0.43 0.43 0.43 0.43 Obj val 314 2272 3971 4292 4294 3PH-A Time (s) 1.1 1.1 1.2 1.2 1.8 Obj val 314 2164 2662 4027 4933 3PH-V Time (s) 1.0 1.2 1.2 1.2 1.1 Obj val 314 2339 2662 4027 4933 Optimal Solution Time (s) 80 1221 2816 1085 6120 Obj val 314 2440 4025 4781 4933

Table 7.7. Percentage Error of Four Heuristics for Five-Airport Problem. Percentage of the total demand 0% 25% 50% 75% 100% Percentage error AGH 0.00% 7.50% 6.09% 7.32% 0.00% VGH 0.00% 6.89% 1.34% 10.23% 12.95% 3PH-A 0.00% 11.31% 33.86% 15.77% 0.00% 3PH-V 0.00% 4.14% 33.86% 15.77% 0.00%

98

AGH 80% 70% 60% Percentage Error 50% 40% 30% 20% 10% 0% 0% 25%

VGH

3PH-A

3PH-V

50%

75%

100%

Number of new devices (% of total demand)

Figure 7.2. Percentage Error of Four Heuristics for the Five-Airport Problem. According to Table 7.5, the computational time of the four methods are relatively small compared to the computational time of the exact algorithm. From Table 7.6, for the five-airport problem, the solution error of AGH is always less than 10% of the optimal solution, giving the best overall performance among four methods. VGH performance is acceptable, but it is not able to obtain the optimal solution in the two extreme cases. The 3PH-A and 3PH-V algorithms obtain the optimal value in the two extreme cases, but their performance declines as the percentage of the total demand is further from the two extreme points (0% and 100%).

99

7.4.2 Large Size Problem The main objective of this research is to study the screening devices allocation problem for both new and existing devices. Therefore, for large size problems, two factors are varied in the experiment: 1) the number of new devices available to be allocated, 2) the number of device types. The large-size problem consists of 100 airports, and three to seven device types. The performances of the heuristics are benchmarked by computing the percentage deviation from upper bound instead of optimal solution, since it was not possible to compute it. The percentage deviation is computed as Percentage Deviation =

UB ? HUE × 100% , UB

where UB is the upper bound of the objective value and HUE is the heuristic solution objective value. The data for 100 airports is randomly generated and is shown Appendix B. The randomly generated data for seven device types is shown in Table 7.8 below.

Table 7.8. Data for Seven Types of Screening Device. Device type 1 2 3 4 5 6 7 Probability of threat detection 0.8 0.7 0.6 0.75 0.95 0.68 0.78 Capacity (bag/device) 150 28 200 40 100 200 120 Operating cost ($/bag) 1 0.83 0.69 0.84 0.94 0.5 0.7 Fixed cost ($/device) 41 12 10 20 60 2 30 Install cost ($/device) 100 40 15 70 250 20 30

100

a) Varying the Number of New Devices Available From the previous experiment, it is suspected that the performance of the heuristics is related to the number of new devices available. Therefore, in this section, we consider the 100-airport problem with various quantities of new devices available. The number of device types is fixed at five types (first five types from Table 7.8). There are two to three classes of passenger for each airport. The number of new devices ranges from 0% to 100% of the total demand of new devices. The total demand is obtained by assuming that there are an unlimited number of new devices and solving each airport individually, as in Section 7.3. The result of the experiment is shown in Table 7.9.

Table 7.9. Comparison of Four Heuristics versus Number of Various New Devices for 100-Airport Problem. Percentage demand 0% 10% 25% 35% 50% 65% 75% 90% 100% (s) 6 10 11 15 17 23 23 25 28 AGH value 33,040 42,075 57,181 64,720 75,615 87,029 95,441 108,340 116,030 (s) 2.3 2.4 2.3 2.4 2.3 2.4 2.4 2.3 2.4 VGH value 32,946 66,948 83,681 93,554 99,685 105,100 105,790 106,240 105,950 (s) 13 29 31 57 36 52 75 38 35 3PH-A value 33,040 52,030 69,934 79,715 93,229 105,210 110,360 115,140 116,030 (s) 15 26 29 54 33 52 73 37 36 3PH-V Objective value 33,040 53,179 71,178 81,342 94,111 105,730 110,310 114,920 116,030

of the total Time Objective Time Objective Time Objective Time

101

According to Table 7.9, the computational times of all heuristics are less than one minute, which is acceptable for a large size problem. The computational time of VGH is stable and also the lowest among the four heuristics. The AGH time increases as the number of new devices increases. For both 3-Phase approaches, the computational times are unpredictable. Considering the two extreme cases (0% and 100%), the AGH, 3PH-A, and 3PH-V obtain the same result, and they are optimal (as they are equal to the upper bound, see Table 7.10), as we expect. Both 3-phase approaches (3PH-A and 3PH-V) outperform the airport greedy heuristic (AGH). These three heuristics select airports to be optimized individually. Therefore, it implies that, for a large number of airports, the rate of improvement (3-Phase selection rule) is a better selection criterion than the number of bags (AGH selection rule). It can be seen from Table 7.9 that when the percentage of the total demand is less than 65%, the VGH performs the best among the four heuristics. For the problem with more than 65% of the total demand, 3PH-A and 3PH-V obtain higher objective values than the other heuristics do. Table 7.10 shows the percentage deviation from the upper bound of the four heuristics. The graph of percentage deviation versus percentage of the total demand is illustrated in Figure 7.3.

102

Table 7.10. Percentage Deviation from Upper Bound versus Percentage of the Total Demand. Percentage of the total demand 0% 10% 25% 35% 50% 65% 75% 90% 100% AGH 0.00% 63.73% 50.70% 44.22% 34.83% 24.99% 17.74% 6.63% 0.00% Percentage deviation from upper bound VGH 0.28% 42.28% 27.85% 19.36% 14.09% 9.42% 8.83% 8.44% 8.69% 3PH-A 0.00% 55.14% 39.71% 31.29% 19.65% 9.33% 4.89% 0.77% 0.00% 3PH-V 0.00% 54.15% 38.63% 29.89% 18.89% 8.88% 4.93% 0.96% 0.00%

The VGH, 3PH-A, and 3PH-V obtain less than 10% of relative deviation from the upper bound for problems with more than 65% of the total demand (see also Figure 7.3). The VGH has the best performance up to 65% of the total demand. Overall percentage deviation from the upper bound is remarkably higher when the percentage of the total demand is between 10%-35%. This is not a result of poor heuristic performance, but because of a weak upper bound. To explain this situation, note that the overall upper bound is the sum of each airport upper bound. The more number of airports the weaker the upper bound because each airport upper bound value does not change when the total number of devices available is more than its new device demand, as a result, the overall upper bound does not decrease as long as the number of devices available is greater than the demand from the airport that has the highest demand.

103

AGH 100% % Deviation from upper bound 90% 80% 70% 60% 50% 40% 30% 20% 10% 0% 0% 20% 40%

VGH

3PH-A

3PH-V

60%

80%

100%

120%

Number of new devices (% of total demand)

Figure 7.3. Percentage Deviation from Upper Bound versus Percentage of the total Demand for 100-Airport Problem. For example, suppose there are five airports to consider, with only one type of new device available. The demands for new devices are 1, 2, 3, 5, and 7 devices respectively (denoted S1, S2,…, S5), and the total demand is 18 devices. Each airport upper bound is computed by solving each airport sub-problem with the number of new devices computed as a percentage of the total number of devices (in this example, 18). When computing the upper bound for total demand equals 18 (or 100%), each airport will attain its maximum upper bound, since the number of devices available for each airport is equal to the total demand, and is more than each airport’s demand (18 > S1, S2,…, S5).

104

Moreover the upper bound, denoted UP(18), is equal to the optimal solution, denoted Z(18), in this case. If the percentage of the total demand decreases by 50%, or the number of devices available equals 9, the upper bound does not decrease. It is still equal to the upper bound of the previous problem, because the number of new devices available is greater than the demand from each airport. However the optimal solution, Z(9), must be less than Z(18) of the previous example, because not all airports will attain their maximum value since there are not enough resources to satisfy all demand. However, UP(9) is equal to UP(18), therefore it is evident that the upper bound, UP(9) is weaker (too high) compared to the optimal solution. Now, if the percentage of the total demand available is 35%, the number of new devices available will be = 18 x 0.35 ≈ 6. The overall upper bound will decrease slightly, as the upper bounds of airports with demand equal to 1, 2, 3, and 5 still attain their maximum value, but the upper bound of the airport with demand equals 7 is decreased (possibly just a small amount). However, the optimal solution should decrease faster than the upper bound does, as now there are at least two airports that could not obtain the maximum value in the optimal solution; one from the airport with demand equals 7 and one or more from the remaining airports which need 1+2+3+5 = 11 and might not receive enough new devices in the optimal solution. In brief, the lower the percentage of the total demand available and the larger the number of airports, the weaker the upper bound. The graph of the upper bound value is depicted in Figure 7.4. It shows that the upper bound is almost stable when the percentage of the total demand ranges from 10%-100%.

105

140,000 120,000 Upper bound 100,000 80,000 60,000 40,000 20,000 0 0% 20% 40% 60% 80% 100% 120% Number of new devices (% of total demand)

Figure 7.4. Upper Bound versus Number of the New Device Available.

b) Varying by the Number of Device Types To compare the performance of the heuristics with various numbers of device types, we create the test problem which consists of 100 airports, two to three classes of passengers for each airport, the number of new devices available equal to 65% of the total demand, and the number of device types is varied from 3 to 7. The percentage of the total demand is chosen to be 65 because, according to previous results, the solution qualities of the VGH, 3PH-A, and 3PH-V are roughly within the same range when the percentage of the total demand is equal to 65. The total demand for three test problems is computed and shown in Table 7.11.

106

Table 7.11. Total Demand for Three Test Problems. Number of device type 3 (first three) 5 (first five) 7 (all seven) Total demand D1 860 542 90 D2 690 76 6 D3 312 218 238 89 7 259 75 993 1242 D4 D5 D6 D7

It is difficult to compare different device types as there are many factors to consider; cost, capacity, or performance. It is interesting to see that the change in the total demand when the number of device types varies may imply the desirability of different device types. From Table 7.11, when there are only three device types to consider the demand for device type one and two are equal to 860 and 690, respectively. When the number of device types is set to seven, the demand for device type six and seven significantly arises while demand for device type one and two decreases dramatically to 90 and 6 devices, respectively. The total demand has shifted from device type one and two to device type six and seven. It implies that device type one and two are less desirable than device type six and seven. Even though device type six and seven have a lower probability of threat detection, they are faster and more cost-effective than device type one and two. Therefore, when we have more choices of device types, it is better to select type six and type seven rather than type one and type two. The experiment is conducted, and the computational time of the four heuristics is shown in Table 7.12. It can be seen that the computational times of the AGH, 3PH-A, and

107

3PH-V noticeably increase as the number of device types increases. To explain this, we know that as the number of device types increases the number of possible solutions exponentially increases, and the number of computational steps of these three heuristics depends on the total number of possible solutions, therefore the computational time increases exponentially. For the VGH, the computational time gradually increases since it is a polynomial time algorithm (the number of computational steps depends on the number of variables).

Table 7.12. Computational Time of Four Heuristic versus the Number of Device Types for 100-Airport Problem. Number of device type 3 5 7 AGH 2.6 24 197 Computational time (second) VGH 1.3 2.3 3.9 3PH-A 9.1 56 243 3PH-V 9.09 52 239

Table 7.13 shows the objective value and the upper bound of each problem. It can be seen that the objective function tends to increase as the number of device types increases. Therefore, we can assume that more choices for devices leads to better security. Note, the VGH outperforms the others when the number of device types equals three. When the number of device types is set to five and seven, VGH, 3PH-A, and 3PHV are in the same range. AGH gives the lowest objective value for all three number of device types. The 3-Phase approaches perform better than VGH when there are three

108

device types, however when there are seven types of devices, VGH obtains better objective function values.

Table 7.13. Objective Value of the Four Heuristics and the Upper Bound versus the Number of Device Types for 100-Airport Problem. Number of device type 3 5 7 AGH 80,495 87,029 92,124 VGH 104,360 104,500 116,580 Objective value 3PH-A 99,821 105,490 113,050 3PH-V 99,788 106,380 115,160 Upper Bound 113,390 116,030 123,800

Consider the percentage deviation from the upper bound in Table 7.14 and Figure 7.5, VGH, 3PH-A, and 3PH-V obtain acceptable objective values, as they are close to the upper bound. We can expect that their objective function values are exceptionally close to the optimal objective function value.

Table 7.14. Percentage Deviation from Upper Bound of Four Heuristics versus the Number of Device Types for 100-Airport Problem. Number of device type 3 5

7

Percentage deviation from upper bound AGH 29.01% 24.99% 25.59% VGH 7.96% 9.94% 5.83% 3PH-A 11.97% 9.08% 8.68% 3PH-V 12.00% 8.32% 6.98%

109

AGH

VGH

3PH-A

3PH-V

50%

Percentage deviation from upper bound

45% 40% 35% 30% 25% 20% 15% 10% 5% 0% 1 3 5 Number of device types 7 9

Figure 7.5. Percentage Deviation from Upper Bound of Four Heuristics versus the Number of Device Types. c) Replicating the Data Set In this section, we experiment thirty sets of randomly generated airport data in order to validate some previous findings and to observe the performance accuracy of the four heuristics. The data consists of 100 airports with two classes of passengers, four types of devices, and the number of available devices is equal to 65% of the total demand. The 2-class scenario is motivated by the available information of the currently employed passenger prescreening system, Secure Flight, and the recommendation from the previous work (McLay et al., 2006). The system divides passengers into two groups; selectee, and ‘non-selectee’ (Singel, 2004). The number of four device types is chosen as currently there are three major types of carry-on baggage screening, and the fourth type

110

will represent the latest technology of screening device. Therefore, for all data sets, the number of existing device type four for all airports will be set to zero. Thirty sets of 100 airports are randomly generated as discussed in Chapter 4. The percentage deviation from the upper bound of each problem by heuristic is shown in Table 7.15 below.

Table 7.15. Percentage Deviation for Thirty Data Set of Four Heuristic with Number of Available Devices = 65% of the Total Demand. Data set 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 AGH 15.51% 14.71% 15.21% 17.18% 11.69% 15.10% 13.70% 12.51% 13.65% 19.03% 12.43% 11.92% 14.24% 11.54% 11.31% 17.47% 16.79% 16.85% 19.25% 15.13% 12.63% 15.89% 13.21% 16.36% 11.73% 14.24% 15.00% 15.57% 14.40% 14.60% VGH 7.18% 6.44% 8.02% 6.42% 7.86% 6.12% 5.75% 5.17% 7.93% 6.78% 6.05% 5.74% 6.85% 6.56% 5.65% 5.43% 6.99% 6.59% 6.12% 6.73% 6.07% 6.85% 7.55% 6.56% 6.62% 5.63% 7.64% 7.66% 6.94% 7.18% 111 3PH-A 6.47% 6.90% 6.46% 7.05% 5.96% 5.99% 6.68% 5.89% 6.79% 6.61% 6.61% 6.49% 7.00% 6.39% 5.87% 6.94% 7.10% 7.01% 7.23% 6.02% 6.31% 6.18% 5.58% 6.55% 6.99% 6.12% 6.78% 6.30% 6.73% 8.08% 3PH-V 6.46% 6.82% 6.46% 7.20% 5.92% 5.87% 7.15% 5.89% 7.17% 6.57% 6.53% 6.49% 7.31% 6.38% 5.83% 6.64% 7.68% 7.12% 7.13% 6.47% 6.30% 6.43% 6.25% 6.70% 7.24% 6.12% 6.72% 6.33% 6.73% 8.08%

Table 7.16. Statistic Data of the Experimental Result in Table 7.15. AGH Average Standard Deviation Maximum Minimum 14.63% 2.15% 19.25% 11.31% VGH 6.64% 0.77% 8.02% 5.17% 3PH-A 6.57% 0.51% 8.08% 5.58% 3PH-V 6.67% 0.54% 8.08% 5.83%

Table 7.16 above presents the statistical data of the experiment. It shows that, on average, the VGH algorithm performs the best among the four Heuristics. The 3PH-A and 3PH-V performances are nearly indistinguishable, and both heuristics’ performances are obviously better than AGH, as expected. In order to continue doing the statistical analysis, a normality test for each heuristic’s replication results must be conducted. The SAS statistical software was used to perform the normality test by using the ShapiroWilk W statistic. The SAS code is given in Appendix C, and tests the following hypothesis: H0 : data are from a normal distribution H1 : data are not from a normal distribution. Table 7.17 shows the P value for each heuristic. The term “P value” refers to the critical level, the descriptive level of significance, or the probability. It indicates the smallest probability level at which we could have preset the level of significant of the test, alpha (α ), and the null hypothesis is still true (Milton and Arnold, 1995). Since all P

112

values are large, we are unable to reject the null hypothesis that the data are drawn from a normal distribution.

Table 7.17. Shapiro-Wilk P Value for Normality Assumption of Four Heuristics. AGH Shapiro-Wilk P value 0.4088 VGH 0.6057 3PH-A 0.4036 3PH-V 0.2369

As we are unable to reject the normality assumption, the next step is to construct a 95% confidence interval for the mean. Since the data is normally distributed and variance is unknown, the T- distribution, Student-t, is used to construct a confidence interval. The confidence interval is calculated as X ± tα / 2 S n , where X , S, n are a sample mean,

standard deviation, and size, respectively, and tα is the point associated with the Tdistribution curve such that the area to the right of the point is equal to α. The 95% confidence intervals of the four heuristics are computed where tα/2 = 2.042, and are given in 7.18 below.

Table 7.18. The 95% Confidence of Four Heuristics’ Percentage Deviation. Lower Bound 13.83% 6.35% 6.38% 6.46% Mean 14.63% 6.64% 6.57% 6.67% Upper Bound 15.43% 6.92% 6.76% 6.87%

AGH VGH 3PH-A 3PH-V

113

From previous experiments, it is suspected that the VGH, 3PH-A, and 3PH-V performances are in the same range when the number of new devices available is equal to 65% of the total demand. A statistical methodology called analysis of variance (ANOVA) was used to test whether the means of these three heuristics are equal. It can be seen that the result of the experiment in Table 7.15 can be considered as a randomized complete block design, where each data set constitutes a block and each heuristic constitutes a treatment. The model of the experiment can be written as

Yij = ? + τ i + β j + ε ij

where

? is an overall mean effect,

τ j is a treatment, or heuristic, effect,

β i is a block, or data set, effect,

and, ε ij is a random deviation attributed to unexplained sources. Before the analysis of variance is conducted, the following assumptions must be validated (Hicks and Turner, 1999): 1. Normal distribution. 2. Independence random samples. 3. Equal population variances. In order to verify the normality assumption, SAS was used to plot the histogram of the residual and perform the normality test by using the Shapiro- Wilk W statistic. The residual is an error term, or what remains after the estimated effects in the model have been subtracted from the values of the response variable. The SAS code is provided in Appendix C, and tests the following hypothesis:

114

H0 : residual is from normal distribution H1 : residual is from normal distribution. Since the P value of the Shapiro- Wilk W statistic, 0.151, is larger than 0.05, and the histogram from SAS software shown in Figure 7.6 is resemble reasonably symmetric and approximately the bell-shape curve, therefore we conclude that the data are from a normal distribution.

The UNIVARIATE Procedure Variable: eresid Stem 9 8 7 6 5 4 3 2 1 0 -0 -1 -2 -3 -4 -5 -6 -7 -8 -9 -10 -11 -12 Leaf 2 3 029 1 166 167789 344467 1255899 0011247778 025 99555331000 764432111 7666653222110 855100 66 88 62 9 83 # 1 1 3 1 3 6 6 7 10 3 11 9 13 6 2 2 2 1 2 Boxplot | | | | | | | +-----+ | | | + | *-----* | | +-----+ | | | | | | 0 0

6 1 ----+----+----+----+ Multiply Stem.Leaf by 10**-3

Figure 7.6. Histogram of Residuals for Single Factor Model.

115

The GLM (General Linear Model) procedure from SAS provides the observed value of the Durbin-Watson (DW) statistic. Since the DW value of 2.4 is large (greater than 1.7), it supports an assumption of independence random sample (Hicks and Turner, 1999). To check the assumption of equal variances, the residuals are plotted against each heuristic. The SAS results from Figure 7.7 do not show the strong evidence to suspect the assumption of equal variance. Therefore, we conclude that variances are equal within treatments. As the model assumptions are satisfied, we are able to compare the three heuristics, VGH, 3PH-A, and 3PH-V. The SAS procedure PROC ANOVA was coded and given in Appendix C. We test the null hypothesis that there are no differences among these heuristics. That is, H0 : all three means are equal H1 : at least one pair of means are not equal. Figure 7.8 shows the selected SAS results of the hypothesis testing. Since the P value of 0.7564 is large, we are unable to reject the null hypothesis and conclude that the three heuristics are not significantly different in this case. From previous experiments, it is suspected that the heuristics 3-PHA and 3-PHV are better than VGH when the percentage of the total demand is large. Therefore, the experiment is repeated with the same 30 sets of airport data, but the percentage of the total demand is changed from 65% to 85%. The Shapiro-Wilk P value, mean, standard deviation, and 95% confidence interval of the four Heuristics are shown in Table 7.19.

116

Plot of eresid*Heu.

Legend: A = 1 obs, B = 2 obs, etc.

eresid ? ? 0.010 ? ?A ? A ?A ?A A ?A ? A A 0.005 ?C A A ? B ?D A A ?B A A ?A A A ?B A D ?A A B 0.000 ? E C ?A B C ?C C B ? D D ?A C D ?A A A ? A -0.005 ?A ?A A ?A A ? ? ?A ? -0.010 ?A ?A ? ? ?A ? ? -0.015 ? ? ????????????????????????????????????????????????????????????????????????? 1 2 3 Heu

Figure 7.7. Residuals versus Heuristic for Single Factor Model.

117

The ANOVA Procedure Dependent Variable: PER Sum of Squares 0.00185153 0.00152230 0.00337383

Source Model Error Corrected Total

DF 31 58 89

Mean Square 0.00005973 0.00002625

F Value 2.28

Pr > F 0.0034

R-Square 0.548792

Coeff Var 0.548656

Root MSE 0.005123

PER Mean 0.933762

Source HEU PROB

DF 2 29

Anova SS 0.00001472 0.00183681

Mean Square 0.00000736 0.00006334

F Value 0.28 2.41

Pr > F 0.7564 0.0022

Figure 7.8. ANOVA Result for Comparing VGH, 3PH-A, and 3PH-V.

Table 7.19. Experimental Result of Four Heuristics when the Percentage of the Total Demand Equals 85. Mean AGH VGH 5.97% 5.66% Standard Deviation Shapiro-Wilk P Value Lower Bound Upper Bound 0.96% 1.05% 0.18% 0.23% 0.4523 0.6571 0.4000 0.2705 5.62% 5.27% 2.28% 2.29% 6.33% 6.05% 2.42% 2.46%

3PH-A 2.35% 3PH-V 2.38%

118

It can be seen that the 3PH-A and 3PH-V performances are very close in terms of mean, standard deviation and the confidence interval. Therefore, we select only one from two of the 3-Phase approaches and compare with the VGH. The algorithm 3PH-V is arbitrary selected. The results of both heuristics are considered as paired data, since one problem set generates two results, one for each heuristic. The difference, D, between each problem set is calculated by Di = 3PH-Vi - VGHi , for i = 1, 2, 3, . . ., 30. To test whether 3PH-V outperforms VGH, the following hypothesis is to be tested: H0 : ?D = 0 H1 : ?D > 0. The SAS code for performing the normality testing of a random sample, D, and a paired T test for above hypothesis is given in Appendix C. Figure 7.9 shows the selected result of the SAS program. It can be seen that, the P value of Shapiro-Wilk test is large, therefore, the difference between two data sets is assumed to be normally distributed. From the MEANS Procedure, the observed value of the test statistic, t value, is equal to 16.83, giving the P value less than 0.001. Since this probability is very small, the null hypothesis is rejected, and we can conclude that, on average, the 3PH-V outperforms the VGH when the 85% of the total demand is met.

119

The UNIVARIATE Procedure Variable: D Moments N Mean Std Deviation Skewness Uncorrected SS Coeff Variation 30 0.0327993 0.01067374 -0.0345609 0.03557776 32.5425796 Sum Weights Sum Observations Variance Kurtosis Corrected SS Std Error Mean 30 0.98397912 0.00011393 -0.8783024 0.00330393 0.00194875

Tests for Normality Test Shapiro-Wilk Kolmogorov-Smirnov Cramer-von Mises Anderson-Darling --Statistic--W D W-Sq A-Sq 0.970867 0.094766 0.035562 0.242397 -----p Value-----Pr Pr Pr Pr < > > > W D W-Sq A-Sq 0.5632 >0.1500 >0.2500 >0.2500

The MEANS Procedure Analysis Variable : D Mean Std Dev Std Error t Value Pr > |t| ??????????????????????????????????????????????????????????????????? 0.0327993 0.0106737 0.0019487 16.83 <.0001 ???????????????????????????????????????????????????????????????????

Figure 7.9. SAS Result for Paired Data.

d) Factorial Experiment In this section, three-factor experiment is conducted in order to study the interaction among the number of device types, number of devices available, and the four heuristics. The percentage deviation from the upper bound of the four heuristics is investigated at three levels of device types (3, 5, and 7) and three levels of the number of devices available (25%, 50%, and 75%). There are twenty four replications (data set) for 120

each cell (combination). The data set consists of 100 airports, and each airport has two classes of passengers. The experiment is considered to be three-factor randomized complete block design, where replication is considered as a block and has no interaction with other factors. The model is also a fixed-effects model, which means the factor levels are pre-determined specifically because they suit the particular interest. There are 3x3x4 = 36 treatment combinations in this experiment. The model of the experiment is written as

Yijkl = ? + β i + H j + Tk + Rl + H j Tk + H j Rl + Tk Rl + H j Tk Rl + ε ijkl

where

? is an overall mean effect,

β i is a block, or data set, effect,

H j is a heuristic effect, T j is a number of device types effect, R j is a number of devices available effect,

and, ε ij is a random deviation attributed to unexplained sources.

Before performing the analysis of variance (ANOVA) the following assumptions must be satisfied: 1. Normal distribution. 2. Independence random samples. 3. Equal population variances.

121

The three-factor experiment was conducted and the percentage deviation was calculated for each combination. SAS was used to plot the histogram, and perform the normality test by using the Kolmogorov-Smirnov test statistic instead of Shapiro- Wilk W statistic as the sample size exceeds 50 (Hicks and Turner, 1999). The SAS code is provided in Appendix C, and tests the following hypothesis: H0 : residual is from normal distribution H1 : residual is from normal distribution. From SAS results, the null hypothesis is rejected in favor of an alternative hypothesis, and we conclude that the original data (percentage deviation) does not conform to the normal distribution. Therefore, the data is transformed by the following function:

(1 ? Y )4.3 , T (Y ) =

4 .3

where Y is the percentage deviation of the original data, and T(Y) is the transformed data. According to SAS results, the P value of the Kolmogorov-Smirnov test statistic for transformed data is larger than 0.05, and the histogram from SAS software shown in Figure 7.10 is reasonably symmetric and shown to have a normal curve. Hence we are unable to reject the null hypothesis and conclude that the data are from a normal distribution.

122

Figure 7.10. Histogram of Residuals for Three-Factor Model.

The Durbin-Watson (DW) test statistic provided by the GLM procedure from SAS is equal to 1.79 which exceeds the threshold value of 1.70. Therefore the assumption of independence random sample is satisfied. Figures 7.11(a) to 7.11(c) separately depict a plot of residuals against each factor. They show no strong evidence against homogeneity of the variances. Therefore, we conclude that variances of residuals are equal within factors, and all the three model assumptions are satisfied.

123

Plot of eresid*Type.

Legend: A = 1 obs, B = 2 obs, etc.

eresid ? 0.03 ? ? ? ? ? ? ?A A 0.02 ?B B ?A B B ?B B B ?A C ?A A ?B A ?D D B 0.01 ?E E B ?J L I ?M J K ?V N T ?S O Z ?Z X Z ?W Z Y 0.00 ?W V Z ?T Y V ?Z V R ?X Q X ?O W R ?P F Q ?I K J -0.01 ?D E H ?H F A ?C E ?C B ?A A B ? A ? B -0.02 ?A ? ? ?A ? ? ? -0.03 ? ????????????????????????????????????????????????????????????????????????? 1 2 3 Type

Figure 7.11. Residuals versus Factor. (a) Number of Device Types.

124

Plot of eresid*Res.

Legend: A = 1 obs, B = 2 obs, etc.

eresid ? 0.03 ? ? ? ? ? ? ? B 0.02 ?B A A ?B A B ?C B A ?B A A ? A A ?B A ?D C C 0.01 ?B C G ?L J I ?H M M ?W P Q ?P U W ?Y Z Y ?Y T Z 0.00 ?Z Z V ?V W V ?R Z W ?Z P R ?S T Q ?P K L ?J K I -0.01 ?D H E ?J B C ?C C B ? A D ? A C ? A ? A A -0.02 ? A ? ? ? A ? ? ? -0.03 ? ????????????????????????????????????????????????????????????????????????? 1 2 3 Res

Figure 7.11. Continued. (b) Number of Devices Available.

125

Plot of eresid*ALGO.

Legend: A = 1 obs, B = 2 obs, etc.

eresid ? 0.03 ? ? ? ? ? ? ? A A 0.02 ? A A B ? B B A ? B A C ? A A B ? B ? A A A ? C C B B 0.01 ? E D A B ? I I E H ? H M E H ? S O M I ? O L N S ? Q R T V ? Q U Z Q 0.00 ? R T V P ? T R O N ? G P V W ? U I S P ? R Q J K ? G I M J ? I H C J -0.01 ? F D A F ? B E D D ? B C B A ? B B A ? B A A ? A ? A A -0.02 ? A ? ? ? A ? ? ? -0.03 ? ???????????????????????????????????????????????????????????? 1 2 3 4 ALGO

Figure 7.11. Continued. (c) Heuristic.

126

As the assumptions of normal distributions, independent random samples, and equal variances are valid, we perform the analysis of variance to study the effect of treatments and interactions. The GLM procedure from SAS was used to investigate this analysis. SAS code is provided in Appendix C, and the result of the experiment is shown in Figure 7.12. It can be seen that all of the factors including the interaction have P value less than 0.0001, therefore it can be concluded that the number of device types, the number of devices available, the heuristic, and interactions significantly have an effect on the model performance.

The GLM Procedure Dependent Variable: Per Sum of Squares 2.18678897 0.03486879 2.22165776

Source Model Error Corrected Total

DF 58 783 841

Mean Square 0.03770326 0.00004453

F Value 846.65

Pr > F <.0001

R-Square 0.984305

Coeff Var 5.392508

Root MSE 0.006673

Per Mean 0.123750

Source Type Res Type*Res ALGO Type*ALGO Res*ALGO Type*Res*ALGO B

DF 2 2 4 3 6 6 12 23

Type I SS 0.01183978 1.64385241 0.00639924 0.41088490 0.01910562 0.04908143 0.00908164 0.03654396

Mean Square 0.00591989 0.82192620 0.00159981 0.13696163 0.00318427 0.00818024 0.00075680 0.00158887

F Value 132.93 18456.9 35.92 3075.56 71.50 183.69 16.99 35.68

Pr > F <.0001 <.0001 <.0001 <.0001 <.0001 <.0001 <.0001 <.0001

Figure 7.12. ANOVA Result for Three-Factor Model from GLM Procedure. 127

It is of interest to identify which heuristic performs the best in each combination of the number of device types and the number of devices available. The StudentNewman-Keuls range test (SNK) is used to compare the four heuristics in each combination. The results are also verified by Duncan’s Multiple Range Test. The example of SAS code for each combination is given in Appendix C. Table 7.20 presents the summary of the results. It can be seen that, overall, the VGH’s performance is significantly better than others heuristics’ performance. Both of the 3-Phase algorithms perform the best when the percentage of total demand equals to 75% and there are 3 and 5 types of device.

Table 7.20. Best Heuristic in Various Conditions. Percentage of total demand 25% Number of device types 3 5 7 VGH VGH VGH 50% VGH VGH VGH 75% 3PH-A & V 3PH-A & V VGH

128

7.5 Chapter Summary In this chapter, experiments were conducted to analyze the performance of the methodologies in terms of solution quality and computational time. The solution quality is based on the percentage error for a small size problem, and percentage deviation from an upper bound of the optimal solution, as the optimal solution cannot be obtained in a reasonable time for a large-size problem. Experimental results show that the performances of the methodologies are related to the number of new devices available and the number of device types to consider. In terms of computational time, the VGH is the best methodology, as expected, since it is the only algorithm that has polynomial computational time. In general, the 3-phase approaches give higher objective function values than AGH does. This means that the rate of improvement is a better selection criterion when one wants to select individual airports to be iteratively optimized. In most cases, the VGH outperforms AGH, 3PH-A, and 3PH-V. However, VGH cannot guarantee the optimal solution in the two extreme cases (no new devices available and an unconstrained number of new devices), while the others always obtain the optimal solution for the extreme cases. When the number of new devices available is insufficient but still greater than approximately 65% of the total demand, the 3PH-A and 3PH-V usually perform better than VGH as well as AGH. In contrast, if the number of new devices available is lower than 65% of the total demand, the VGH outperforms other methodologies. The overall quality of all four heuristics is improved as the number of device types increases.

129

CHAPTER VIII CONCLUSIONS AND FUTURE RESEARCH

8.1 Conclusions In the preceding chapters, we presented the development of an explosive screening device allocation model for airport security systems. In Chapter 1, we introduced the background of aviation security and details of passenger screening. Explosive screening devices for passenger screening are grouped into three categories: passengers, carry-on baggage, and checked baggage. We addressed existing literature related to this research in Chapter 2. The literature survey was divided into three related areas: aviation security, complexity theory, and optimization problems. The optimization model for single airport (station) security systems is developed in Chapter 3, where we focused on carry-on baggage screening. Before boarding an aircraft, passengers are grouped into several classes according to their level of risk. The model goal is to maximize the probability of threat detection. The model also accounts for capacity and budget constraints. Complexity analysis shows that the explosive screening device allocation model is NP-hard by presenting a polynomial transformation from a known NP-complete problem; the Binary Knapsack Problem. In Chapter 4, we discuss Lawler and Belle’s partial enumeration algorithm (Lawler and Belle, 1966) which is the primary method used to obtain the optimal solution of the single-airport problem. The modification of the model to be consistent with Lawler and Belle’s standard form is also shown. We provide numerical examples and the results of the algorithm are shown.

130

A total enumeration algorithm is used to obtain results for several small problems to verify the accuracy and the speed of the algorithm. In Chapter 5, to address the multiple airports problem, the model is expanded from a single-airport model to include multiple airports in one optimization problem. The primary objective of multiple-airport model is to maximize the total security level of all airports under a set of realistic constraints. The objective function is derived from the probability of threat detection, passenger average assessed threat value, and number of bags. The constraints account for the total budget, system throughput, and new device resource availability. The multiple-airport model is a large scale optimization problem with nonlinear objective function and multiple constraints. The model is considered to be NP-hard problem, as the single-airport, restricted version, is shown to be NP-hard problem. For that reason, an optimal solution for a multiple-airport problem may not be obtainable in a reasonable time. Moreover, most of the data in the model are estimated. Solving the model to optimality with approximate data may not be worthwhile in this situation. For that reason, a simple and less time-consuming heuristic method is worthwhile for obtaining the near optimal solution instead of a global optimal solution. Four heuristics are presented in Chapter 6. They are 1) Airport Greedy Heuristic (AGH), 2) Variable Greedy Heuristic (VGH), 3) Three-Phase Approach Combined with Airport Greedy Heuristic (3PH-A), and 4) Three-Phase Approach Combined with Variable Greedy Heuristic (3PH-V). Experiments are conducted in Chapter 7 to examine the heuristics’ solution quality and computational time. For small size problems, the heuristics are compared to

131

the optimal solution. For large size problems, the optimal solution is not obtainable in a reasonable time, so the upper bound is used instead of the optimal solution for large size problems. It is shown in small size problems that the AGH, 3PH-A, and 3PH-V obtain the optimal solution in two extreme cases (no new device available and sufficient new device available). The results of the experiment also indicate that the VGH is effective for most problems in term of both solution quality and computational time. However, when the number of new devices available is high, i.e. more than 70% of total demand, and the number of device types equals to 3 and 5, both of the 3-Phase approaches achieve the best solution quality among four heuristics.

8.2 Directions of Future Research Future research can be done in several directions depending on the availability of data, and the objective of the study. Three main directions are discussed in this section. 1. Expand model to include passenger and checked baggage screening: The objective function in the model could be adapted to incorporate two other ways of screening a threat. We know that there are three means that a terrorist can bring a threat onto an aircraft: checked baggage, carry-on baggage, and on their person. By using a weighted average, the probability of threat detection in the objective function may be modified from probability of threat detection for carry-on baggage to all three means of threat carrying. Initially, each way of carrying a threat may be treated equally. However, if we can determine the probability of

132

terrorist means of threat carrying, we can distribute more accurate proportion to each way. 2. False Alarm: Usually, a false alarm is not desirable. It can lead to passenger inconvenience. We might want to control the false-alarm rate in each airport, or possibly want to reduce it. The first option could be accomplished by adding a false alarm constraint for each airport to limit the false-alarm rate to be lower than some predefined value. This restriction will certainly reduce the feasible domain of the solution or lead to problem infeasibility. For the latter option, we can include a false alarm in the objective function or completely replace the probability of threat detection (true alarm) with the probability of false alarm. 3. Applications to other problems: To better conform to real world scenarios, the methodologies presented in this dissertation might be experimented with using real data. The data may consist of approximately 500 airports in one optimization problem. The methodologies may also be applied to other types of optimization problems that have similar structure. The problem might consist of one or more of these characteristic: a block structure, a nonlinear objective function, or binary variables, etc. Examples are reliability, transportation, or many types of allocation problems.

133

REFERENCES

Aggarwal, K.K., Gupta, J.S., and Misra, K.B. (1975). A new heuristic criterion for solving a redundancy optimization. IEEE Transactions on Reliability 24(1), 8687. Barnett, A. (2004). CAPPS II: The foundation of aviation security? Risk Analysis 24, 909-916. Bretthauer, K.M., and Shetty, B. (1995). The nonlinear resource allocation problem. Operations Research 43, 670-683. Bretthauer, K.M., and Shetty, B. (2002). A pegging algorithm for nonlinear resource allocation problem. Computers and Operations Research 29, 505-527. Butler, V., and Poole Jr., R.W. (2002). Rethinking Checked-Baggage Screening. Policy Study 297. Reason Public Policy Institute, Los Angeles, CA. Chern, M.S., and Jan, R.H. (1986). Reliability Optimization Problems with Multiple Constraints. IEEE Transactions on Reliability 35(4), 431-436. Cobb, R.G., and Primo, D.M. (2003). The Plane Truth: Airline Crashes, the Media, and Transportation Policy, Brookings Institution Press, Washington, D.C. Cook, S.A. (1971). The complexity of theorem-Proving procedures. Proceedings of the 3rd Annual ACM Symposium on Theory of Computing, Association for Computing Machinery, New York, 151-158. Floudas, C.A. (1995). Nonlinear and Mixed-Integer Optimization, Oxford University Press, Oxford. Gary, M.R., and Johnson, D.S. (1979). Computers and Intractability: A Guide to the Theory of NP-Completeness. W.H. Freeman and Company, New York, NY. Hall, L. (2000). Encyclopedia of Operations Research and Management Science. 2nd Edition, Springer, New York, NY. Hicks, C.R., and Turner, Jr., K.V. (1999) Fundamental Concepts in the Design of Experiments. 5th Edition, Oxford University Press, New York, NY.

134

Jacobson, S.H., Virta, J.E., Bowman, J.M., Kobza, J.E., and Nestor, J.J. (2003) Modeling aviation baggage screening security systems: a case study. IIE Transactions 35, 259-269. Kellerer, H., Ulrich, P., and Pisinger, D. (2004). Knapsack Problems. Springer, Heidelberg, Germany. Kobza, J.E., and Jacobson, S.H. (1997). Probability models for access security system architectures. Journal of Operational Research Society 48, 255-263. Kovalyov, M.Y. (1996). A rounding technique to construct approximation algorithms for knapsack and partition-type problems. Applied Mathematics and Computer Science 6, 789-801. Lawler, E.L., and Belle, M.D. (1966). A method for solving discrete optimization problems. Operations Research 14, 1098-1112. Lazar Babu, V.L. (2003). Airport Security System Design. Buffalo: Univ. at Buffalo (SUNY). Thesis. Lueker, G.S. (1975). Two NP-complete problems in nonnegative integer programming. Technical Report 178, Computer Science Laboratory, Princeton University, Princeton, NJ. Lusena, C. (1999). Easy versions of the knapsack problem. N.p.;n.p., author URL http://www.users.csbsju.edu/~clusena/ Mao, J.C.T., and Wallingford, B.A. (1968). An extension of Lawler and Belle’s method of discrete optimization with examples from capital budgeting. Management Science 15, B51-B60. Mathur, K., Salkin, H.M., and Mohanty, B.B. (1986). A note on a general non-linear knapsack problem.” Operation Research Letters 5, 79-81. McLay, L.A., Jacobson, S.H., and Kobza, J.E. (2006). Multilevel Passenger Screening Strategies For Aviation Security Systems. Naval Research Logistics (Accepted). Mead, K.M. (2002). Key challenges facing the transportation security administration. CC-2002-180. Washington, DC: U.S. Department of Transportation, Office of the Inspector General. Milton J.S., and Arnold J.C. (1995). Introduction to Probability and Statistics: Principles and Applications for Engineering and the Computing sciences. 3rd Edition. McGraw-Hill, New York, NY.

135

Misra, K.B. (1971) A method of solving redundancy optimization problems. IEEE Transactions on Reliability 20, 117-120. Misra, K.B., and Sharman, J. (1973). Reliability optimization of a system by zero-one programming. Micro-Electronics and Reliability 12, 229-233. Murty, K.G., and Kabadi, S.N. (1987). Some NP-complete problems in quadratic and nonlinear programming. Mathematical Programming 39, 117-129. National Research Council. (1996). Airline Passenger Security Screening: New Technologies and Implementation Issues. National Academy Press Publication NMAB-482-1, Washington, D.C. Nemhauser, G.L., and Wolsey, L.A. (1999) Integer and Combinatorial Optimization. John Wiley & Sons, Inc., New York, NY. Pardolos, P.M. (1996). Continuous approaches to discrete optimization problems In Nonlinear and Optimization and Applications, G. Di Pillo & F. Giannessi, Ed., Plenum Press, New York, NY 313-328. Poole Jr., R.W., and Passantino G. (2003). A risk-based airport security policy. Policy Study 308. Reason Public Policy Institute, Los Angeles, CA. Rausch, J. (2002). Optimizing the deployment schedule of a constrained resource. INFORMS Annual Meeting, Northrop Grumman Information Technology. Ritchie, L. (2004). A Cost-Benefit Analysis of Alternative Device Configurations for Aviation Checked Baggage Security Screening. Lubbock: Texas Tech Univ. Thesis. Royal Air Force Museum. (2002). World Aviation in 1909. [WWW document]. URL http://www.rafmuseum.org.uk/milestones-of-flight/world/1909.html Sharma, J., and Venkateswaran, K.V. (1971). A direct method for maximizing the system reliability. IEEE Transactions on Reliability 28(3), 256-259 (1971). Sharma, U., and Misra, K.B. (1990). An efficient algorithm to solve integer-programming problems in reliability optimization. International Journal of Quality & Reliability Management 7(5), 44-56. Singel, R. (July 14, 2004). Airline Screening System Delayed. Wired News [Online newspaper]. URL http://www.wired.com/news/privacy/ 0,1848,64201,00.html.

136

Tsurkov, V. (2001). Large-scale Optimization- Problems and Methods. Kluwer Academic Publishers, Norwell, MA. U.S. Centennial of Flight Commission.(2004). Aviation Security [WWW document]. URL http://www.centennialofflight.gov/essay/Government_Role/ security/POL18.htm U.S. Aviation Subcommittee. (2001). Aviation Security and the Future of the Aviation Industry (Hearing, 107th Congress, 1st Sess.). Washington, DC: U.S. Government Printing Office. U.S. General Accounting Office. (2003). Airport passenger screening preliminary observations on progress made and challenges remaining. GAO-03-1173. Washington, DC: U.S. Government Printing Office. U.S. Department of Homeland Security (2004) Privacy Act of 1974: System of Records; Secure Flight Test Records. Document TSA-2004-19160-584. Virta J.E., Jacobson S.H., and Kobza J.E. (2002). Outgoing Selectee Rates at Hub Airports. Reliability Engineering and System Safety 76, 155-165.

Virta J.E., Jacobson S.H., and Kobza J.E. (2003). Analyzing the Cost of Screening Selectee and Non-Selectee Baggage. Risk Analysis 23, 897-908.

137

APPENDIX A MATLAB SOURCE CODE

A.1 Single Airport Problem A.1.1 Main Algorithm

% Author: Jackrapong Attagara % Industrial Engineering % Texas Tech University % This code calls Lawler and Bell's algorithm and total enumeration to find % the optimal solution for a single airport problem. clear all; Initialize; tic; LawlerBell; toc; tic; SearchAll; toc; ; function Initialize() % load data from excel and assign value T = xlsread('data.xls'); global G; G = T(1,1); global M; M = T(2,1); global B; B = T(3,1); global N; N = T(5,1:G); global AT; AT = T(6, 1:G); global P; P = T(8, 1:M); global CP; CP = T(9, 1:M); global OC; OC = T(10, 1:M); global FC; FC = T(11, 1:M); global IC; IC = T(12, 1:M); global E; E = T(13, 1:M); global R; R = CalR global X; X = zeros(1, M*G);

139

;

function LawlerBell() disp('=================== Lawler and Belle for single problem==================='); global G M B N AT P CP OC FC IC E R X; bestObj = 1000; bestG11_X =20; bestX = zeros(1, length(X)); overflow = false; overflow = checkOverflow(X); loop_count = 0; skip_count = 0; XStar_Overflow = false; XPlusOne_Overflow = false; while overflow==false neg_X = 1-X; Y = CalY(neg_X, N, CP); S = CalS(Y, E); Temp_Cost = CalCost(neg_X, Y, S, OC, FC, IC, N); G11_X = B - Temp_Cost; newObj = -calculateObj(neg_X, P, R); XStarMinusOne = CalXStarMinusOne(X); neg_X = 1-XStarMinusOne; Y = CalY(neg_X, N, CP); S = CalS(Y, E); Temp_Cost = CalCost(neg_X, Y, S, OC, FC, IC, N); G11_XStarMinusOne = B - Temp_Cost; XStar_Overflow = checkOverflow(XStarMinusOne); XplusOne_Overflow = checkOverflow(X); XStar = function_XplusOne(XStarMinusOne); loop_count = loop_count+1; if newObj >= bestObj %('Rule 1, x = x*') X = XStar; skip_count = skip_count+1; if XStar_Overflow == true overflow = true; end else if G11_XStarMinusOne <= 0 %('Rule 2, x = x*') X = XStar; skip_count = skip_count+1; if XStar_Overflow == true overflow = true; %('overflow 1') end else if G11_X >= 0 %('Rule 3, x=x*')

140

bestX = X; bestObj = newObj; bestG11_X = G11_X; X = XStar; skip_count = skip_count+1; if XStar_Overflow == true overflow = true; %('overflow 2') end else if XPlusOne_Overflow == true overflow = true; %('overflow 3') end %('no rule x=x+1'); X = function_XplusOne(X); end end end if overflow==true disp('overflow, stop'); end XStar_Overflow = false; XplusOne_Overflow = false; end mat_x = CalMatX(1-bestX, length(N), length(CP)); Y = CalY(1-bestX, N, CP); S = CalS(Y, E); bestG11_X = CalCost(1-bestX, Y, S, OC, FC, IC, N); Best = ['' int2str(mat_x) ''] Objective= ['' num2str(-bestObj) ' Cost = ' num2str(bestG11_X) ' '] disp('=================== END LAWLER AND BELLE===================');

function SearchAll() % this function search all possible solution and keep the best solution disp('=================== Search All==================='); global G M B N AT P CP OC FC IC E R X; Y = CalY(X, N, CP); S = CalS(Y, E); Cost = CalCost(X, Y, S, OC, FC, IC, N); Obj = calculateObj(X, P, R); if Cost < B bestX = X; end overflow = false; overflow = checkOverflow(X); while overflow==false X = function_XplusOne(X); newObj = calculateObj(X, P, R); Y = CalY(X, N, CP); S = CalS(Y, E); mat_x = CalMatX(X, length(N), length(CP));

141

newCost = CalCost(X, Y, S, OC, FC, IC, N); if (newObj > Obj) & (newCost <= B) bestX = X; Obj = newObj; Cost = newCost; end overflow = checkOverflow(X); if overflow==true disp('overflow, stop') end end mat_x = CalMatX(bestX, length(N), length(CP)); Best = ['' int2str(mat_x) ''] Objective= [' ' num2str(Obj) ' Cost = ' num2str(Cost) ' '] disp('=================== END ===================');

A.1.2 Public Function

function y = CalY(x, n, cp) % calculate number of devices needed mat_X = CalMatX(x, size(n,2), size(cp,2)); M = zeros(1, length(cp)); for i=1:length(M); M(i) = ceil(n*mat_X(:,i)/cp(i)); end y = M; function y = CalS(n,e) %calculate number of new devices needed L = length(n); M = zeros(1,L); for i = 1:L if n(i)>e(i) %if new device needed M(i)=n(i)-e(i); else % if no new device M(i)=0; end end y = M; function y=checkOverflow(x) % This function check whether the vector is overflow if x == ones(1, length(x)) overflow = true; else overflow = false; end y = overflow; function cost = CalCost(x, y, s, oc, fc, ic, n)

142

%calculate fix + installation + operating cost y_cost = y*fc' + s*ic' ; x_cost =zeros(1,length(fc)); mat_X = CalMatX(x, size(n,2), size(fc,2)); for i=1:length(fc) x_cost(i) = n*mat_X(:,i)*oc(i); end cost = y_cost + sum(x_cost); function y=calculateObj(x, p, r) %calculate the objective function if size(p) ~= size(r) disp('Error in dimension of AT and N !!!'); else m = length(p); %number of machine g = length(r); %number of passenger group L=zeros(1,g); for i=1:1:g; gx = x(1+(i-1)*m : i*m); L(i)=(1-prod(1-gx.*p))*r(i); end end y = sum(L); function y=CalXStarMinusOne(x) %find the binary vector x star minus one form vector x tempx = x; if x == zeros(1,length(x)) tempx = x; else x(end)=x(end)-1; for i=length(x):-1:1 if x(i)< 0 x(i)=1; if i>1 x(i-1)=x(i-1)-1; end end end end y = zeros(1,length(x)); for j=1:length(x) y(j)= tempx(j)|x(j); end ; function y=checkOverflow(x) % This function check whether the vector is overflow (all 1’s) if x == ones(1, length(x)) overflow = true; else overflow = false; end

143

y = overflow;

function y=function_XplusOne(x) %find the vector X plus one from vector x x(end)=x(end)+1; for i=length(x):-1:1 if x(i)>1 x(i)=0; if i>1 x(i-1)=x(i-1)+1; else %disp('overflow from plus one function'); end end end y=x; function y = CalMatX(x,r,c) %change the binary vector to matrix form to display the solution L = length(x); i=1; while i <= L for j=1:r for k=1:c M(j,k)=x(i); i = i+1; end end end y = M;

A.2 Multiple Airport Problem A.2.1 The Airport Greedy Heuristic (AGH) Main Algorithm

tic clear global Airport; KReadData; %read data NumberOfBag = CalNumberOfBag; %calculate total number of bags [xs, idx] = sort(NumberOfBag, 'descend'); global U; global K; global Airport; Total_Obj =0; uRecent = U; for i=1:K j = idx(i); HLawerAndBell(j, uRecent); % same as LawlerBell except the number of available

144

devices, U uRecent = uRecent - Airport(j).bestS; Total_Obj = Total_Obj + Airport(j).H1Obj; i=i+1; end Total_Obj = Total_Obj disp('xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx'); toc

A.2.2 The Variable Greedy Heuristic (VGH) Main Algorithm

tic clear all; KReadData; %read data global Airport; global K; % of airport global P; % probability of threat detection global D; %max # device global CP; % capacity of devices global FC; % device fixed cost global IC; % device install cost global OC;% device operating cost global U; % # of available device maxAirport =K; maxDevice = D; global G; % max # of group Factor = zeros(maxAirport, 3, maxDevice);%prob increase if x=1 Y = zeros(maxAirport, 3, maxDevice);%# of device use if x=1 S = zeros(maxAirport, 3, maxDevice);%# of new device needed if x=1 Cost = zeros(maxAirport, 3, maxDevice);%cost use if x=1 Profit = zeros(maxAirport, 3, maxDevice);%obj increase if x=1 Cost_Ratio = zeros(maxAirport, 3, maxDevice); Pass_So = ones(maxAirport, 3, maxDevice); for i=1:maxAirport for j=1:Airport(i).numGroup for k=1:maxDevice Profit(i,j,k) =P(k)*Ai