Thomas Morris

Wei Gao

INDUSTRIAL CONTROL SYSTEM TRAFFIC DATA SETS FOR INTRUSION DETECTION RESEARCH

Keywords: Industrial control systems, SCADA, intrusion detection, MODBUS

Supervisory control and data acquisition (SCADA) systems monitor and control physical processes associated with the critical infrastructure. Weaknesses in the application layer protocols, however, leave SCADA networks vulnerable to attack. In response, cyber security researchers have developed myriad intrusion detection systems. Researchers primarily rely on unique threat models and the corresponding network traffic data sets to train and validate their intrusion detection systems. This leads to a situation in which researchers cannot independently verify the results, cannot compare the effectiveness of different intrusion detection systems, and cannot adequately validate the ability of intrusion detection systems to detect various classes of attacks. Indeed, a common data set is needed that can be used by researchers to compare intrusion detection approaches and implementations. This paper describes four data sets, which include network traffic, process control and process measurement features from a set of 28 attacks against two laboratoryscale industrial control systems that use the MODBUS application layer protocol. The data sets, which are freely available, enable effective comparisons of intrusion detection solutions for SCADA systems.

Introduction

Supervisory control and data acquisition (SCADA) systems are computerbased process control systems that control and monitor remote physical processes. SCADA systems are strategically important because they are widely used in the critical infrastructure. Several incidents and cyber attacks affecting SCADA systems have been documented; these clearly illustrate the vulnerability of critical infrastructure assets. The reported incidents demonstrate that cyber attacks against SCADA systems can have severe financial impact

System

Threat Model Protocol SRI Modbus [START_REF] Cheung | Using model-based intrusion detection for SCADA networks[END_REF] Access, reconnaissance and attack MODBUS NNIDSCI [START_REF] Linda | Neural network based intrusion detection system for critical infrastructures[END_REF] Traffic from Nmap, Nessus, Metasploit -AKKR-SPRT [START_REF] Yang | Anomaly-based intrusion detection for SCADA systems[END_REF] DoS attacks simulated by Sun servers SNMP IDAEM [START_REF] Oman | Intrusion detection and event monitoring in SCADA networks[END_REF] RTU attacks -Multidimensional CSA [START_REF] Carcano | A multidimensional critical state analysis for detecting intrusions in SCADA systems[END_REF] Simulated attacks on critical states MODBUS SGDIDS [START_REF] Zhang | Distributed intrusion detection system in a multi-layer network architecture of smart grids[END_REF] KDD 99 Cup Data Set -Pattern Detection [START_REF] Valdes | Communication pattern anomaly detection in process control systems[END_REF] Reconnaissance MODBUS KSSM [START_REF] Linda | Improving control system cyberstate awareness using known secure sensor measurements[END_REF] False data injection -Statistical Estimation [START_REF] Rrushi | Detecting anomalies in process control networks[END_REF] Overflow exploits MODBUS RAIM [START_REF] Ten | Anomaly detection for cybersecurity of substations[END_REF] File system and status modification C37.118

and can result in damage that is harmful to humans and the environment.

In 2000, a disgruntled engineer compromised a sewage control system in Maroochy Shire, Australia, causing approximately 264,000 gallons of raw sewage to leak into a nearby river [START_REF] Slay | Lessons learned from the Maroochy water breach[END_REF]. In 2003, the Slammer worm caused a safety monitoring system at the Davis-Besse nuclear plant in Oak Harbor, Ohio to go offline for approximately five hours [START_REF] Poulsen | Slammer worm crashed Ohio nuke plant network, Security-Focus, Symantec[END_REF]. The insidious Stuxnet worm [START_REF] Falliere | Stuxnet Dossier, Version 1.4[END_REF], which was discovered in 2010, targeted nuclear centrifuge system controllers, modifying system behavior by distorting monitored process information and altering control actions. Cyber security researchers have developed numerous intrusion detection systems to detect attacks against SCADA systems. Much of the research uses training and validation data sets created by the same researchers who developed the intrusion detection systems. Indeed, no standardized data set containing normal SCADA network traffic and attack traffic is currently available to researchers. In order to evaluate the performance of data mining and machine learning algorithms for SCADA intrusion detection systems, a network data set used for benchmarking intrusion detection system performance is sorely needed. This paper describes four data sets, which include network traffic, process control and process measurement features from a set of 28 attacks against two laboratory-scale industrial control systems that use the MODBUS application layer protocol. The data sets, which are freely available, enable effective comparisons of intrusion detection solutions for SCADA systems.

Related Work

Several SCADA security researchers have developed intrusion detection systems that monitor network traffic and detect attacks against SCADA systems. Table 1 lists example intrusion detection systems, the threat models they use and the network protocols they analyze. Note that each intrusion detection system uses a unique threat model. Some threat models are based on attacks executed against SCADA laboratory testbeds while others are based on ma-nipulated data sets drawn from other domains. The network protocols also differ; MODBUS is the most common protocol (used in three systems) while the IEEE C37.118 protocol is used in just one system. The remaining systems use threat models with attacks implemented at different network layers.

A noticeable drawback of the research identified in Table 1 is that the threat models only include subsets of attack classes. Not surprisingly, exploit coverage is limited for each of the data sets. Only a few of the threat models consider reconnaissance attacks while some models only include response injection attacks. Indeed, the malicious behavior captured in the data sets is neither consistent nor comprehensive in terms of normal operations and attacks. For this reason, it is difficult to judge the effectiveness of an intrusion detection system against sophisticated attacks. This also leads to a situation in which researchers cannot independently verify intrusion detection results and cannot compare the performance of intrusion detection systems.

Test Bed Description

The data sets described in this paper were captured using a network data logger, which monitored and stored MODBUS traffic from a RS-232 connection. Two laboratory-scale SCADA systems were used: a gas pipeline and water storage tank.

Figure 1 shows the gas pipeline and water storage tank systems along with the associated human machine interfaces (HMIs). The gas pipeline system includes a small airtight pipeline connected to a compressor, a pressure meter and a solenoid-controlled relief valve. The pipeline system attempts to maintain the air pressure in the pipeline using a proportional integral derivative (PID) control scheme.

The water storage tank system includes a tank that holds approximately two liters of water, a manually-operated relief valve to deplete water from the tank, a pump to add water to the tank from an external water source and a meter to measure the water level as percentage of tank capacity. The water storage tank uses an on/off control scheme to maintain the water level between the high (H) and low (L) setpoints. The water storage tank activates an alarm when the water level is above the high alarm setpoint (HH) or below the low alarm setpoint (LL). Detailed descriptions of the functionality of the two systems and their respective components are provided in a separate paper [START_REF] Morris | A control system testbed to validate critical infrastructure protection concepts[END_REF].

A bump-in-the-wire approach was used to capture data logs and to inject attacks. The device was implemented via a C program running on a VMware virtual machine. The virtual machine included two RS-232 serial ports connected to a USB-to-serial converter. The C program monitored each serial port for traffic. Detected traffic was timestamped and recorded in a log file. To facilitate attacks, the C program incorporated hooks to inject, delay, drop and alter network traffic.

Description of Attacks

The data sets presented in this paper include network traffic, process control and process measurement features from normal operations and attacks against the two SCADA systems. The attacks are grouped into four classes: (i) reconnaissance; (ii) response injection; (iii) command injection; and (iv) denial-of-service (DoS).

Reconnaissance Attacks

Reconnaissance attacks gather SCADA system information, map the network architecture and identify device characteristics (e.g., manufacturer, model number, supported network protocols, device address and device memory map). The reconnaissance class of attacks in the data set includes four attacks against MODBUS servers: address scan, function code scan, device identification attack and points scan. The address scan discovers SCADA servers connected to a network by polling for responses from different MODBUS addresses. The function code scan identifies supported MODBUS function codes that can be used by an identified server. The device identification attack allows an attacker to obtain device vendor information, product code and major and minor firmware revisions. The points scan allows the attacker to build a memory map of MODBUS coils, discrete inputs, holding registers and input registers.

Response Injection Attacks

SCADA systems commonly use polling techniques to continuously monitor the state of a remote process. Polling takes the form of a query transmitted from the client to the server followed by a response packet transmitted from the server to the client. State information is provided to a human machine interface for monitoring the process, storing process measurements in a data historian and providing feedback to control loops that measure process parameters and take the appropriate control actions based on the process state. Response injection attacks alter responses from the server to client, providing false system state information.

Response injection attacks are divided into naive malicious response injection (NMRI) attacks and complex malicious response injection (CMRI) attacks. NMRI attacks leverage the ability to inject or alter response packets in a network; however, they lack the ability to obtain information about the underlying process being monitored and controlled. Eight NRMI attacks were used in creating the data sets described in this paper. The naive read payload size attack returns a malicious response with the correct payload size but sets the payload to all zeros, ones or random bits. The invalid read payload size attack returns a malicious response with a length that does not conform to the requested length. The invalid exception code attack returns false error responses to the client after a read command. The negative sensor measurements attack injects negative process measurements; this is problematic because many systems use floating point numbers to represent values that can only be positive. The sensor measurements grossly out-of-bounds attack injects process measurements that are significantly outside the bounds of alarm setpoints. The sporadic sensor measurement injection attack sends false process measurements outside the bounds of the H and L control setpoints while staying within the alarm setpoint range specified by HH and LL. The random sensor measurement injection attack sends random process measurements of gas pipeline pressure or water tank water level.

CMRI attacks attempt to mask the actual state of the physical process and negatively affect feedback control loops. They are more sophisticated than NMRI attacks because they require an in-depth understanding of the targeted system. As such, CMRI attacks are designed to appear like normal process functionality. These attacks can be used to mask alterations to process state perpetrated by malicious command injection attacks. CMRI attacks are more difficult to detect because they project a state of normalcy.

Five CMRI attacks were used to create the data sets. The constant sensor measurement injection attack repeatedly sends malicious packets containing the same measurement to mask the real state of the system. The calculated sensor measurement injection attack sends pre-calculated process measurements. The high frequency measurement injection attack increases the rate of change of a process measurement beyond its normal range. The low frequency measurement injection attack decreases the rate of change of a process measurement below its normal range. A replayed measurement injection attack resends process measurements that were previously sent from the server to a client.

Command Injection Attacks

Command injection attacks inject false control and configuration commands to alter system behavior. The potential impacts of malicious command injections include loss of process control, interruption of device communications, unauthorized modification of device configurations and unauthorized modification of process setpoints. Command injection attacks are divided into malicious state command injection (MSCI) attacks, malicious parameter command injection (MPCI) attacks and malicious function code command injection (MFCI) attacks. Comprehensive descriptions of these attacks are provided in [START_REF] Gao | Cyber Threats, Attacks and Intrusion Detection in Supervisory Control and Data Acquisition Networks[END_REF].

MSCI attacks change the state of the process control system to drive the system from a safe state to a critical state by sending malicious commands to remote field devices. MSCI attacks may involve a single injected command or multiple injected commands. Three MSCI attacks were used to create the data sets. The altered system control scheme attack changes the control mode from automatic to manual and then turns on the compressor or pump to increase the pressure in the pipeline or raise the water level in the water storage tank, respectively. The altered actuator state attack changes the state of an actuator in a system. In the case of the gas pipeline system, this attack includes command injections that turn the compressor on or off, and those that open or close the relief valve; in the case of the water storage tank system, the altered actuator state attack turns the pump on or off. The continuous altered actuator state attack repeatedly changes the actuator states in a system. For example, command packets could be continually transmitted to switch the state of the compressor and pump in the pipeline and storage tank systems, respectively. Additionally, a continuous altered actuator state attack may be used to repeatedly transmit MODBUS write register commands to invert the state of the solenoid that controls the relief valve in the gas pipeline system. MPCI attacks alter programmable logic controller (PLC) field device setpoints. The data sets include two MPCI attacks. The altered control setpoint attack changes the H and L setpoints for the water storage tank while disabling the liquid level alarms. A proportional integral derivative (PID) controller is commonly used in SCADA systems to maintain a desired setpoint by calculating and adjusting for system error; the altered proportional integral derivative parameter attack changes the PID parameters used in the gas pipeline system.

MFCI attacks use built-in protocol functions in a manner different from what was intended. The data sets include four MFCI attacks. The force listen only mode attack causes a MODBUS server to stop transmitting on the network. The restart communications attack sends a command that causes the MODBUS server to restart, leading to a temporary loss of communications. The clear communications event log attack erases the communications event log of the MODBUS server. Finally, the change ASCII input delimiter attack changes the delimiter used for MODBUS ASCII devices.

Denial-of-Service Attacks

Denial-of-service attacks target communications links and system programs in an attempt to exhaust resources. The data sets include two denial-of-service attacks. The invalid cyclic redundancy code (CRC) attack injects a large number of MODBUS packets with incorrect CRC values into a network. The MODBUS master traffic jamming attack uses a non-addressed slave address to continually transmit random data to random destination addresses.

SCADA Traffic and Payload Data Sets

The KDD Cup 1999 Data Set [START_REF] Hettich | The UCI KDD Archive[END_REF] was developed to train and validate intrusion detection systems associated with traditional information technology systems. The use of this common data set by numerous researchers facilitated the independent validation of research results and the comparison of many intrusion detection system approaches. In the area of SCADA security, however, researchers develop their own data sets to test intrusion detection systems because there is a lack of availability and access to SCADA network traffic. Indeed, no standard data set is available that includes normal and attack traffic for a SCADA network that can serve as a benchmark to evaluate and compare SCADA intrusion detection system performance. This section describes a data set that is intended to provide researchers with a common platform to evaluate the performance of data mining and machine learning algorithms designed for SCADA intrusion detection systems. The data set includes different classes of attacks that cover a variety of SCADA system attack scenarios.

The common data set described in this paper has three primary benefits. First, not all researchers have access to SCADA equipment to generate their own data sets; a common data set would enable more researchers to work in the area of SCADA security. Second, a common data set would allow researchers to independently validate the results of other researchers. Third, a common data set would enable the comparison of the performance of different algorithms, leading to better intrusion detection systems.

Data Set Organization

The data sets created as a result of this research effort are stored in the Attribute Relationship File Format (ARFF) for use with the WEKA software [START_REF] Hall | The WEKA data mining software: An update[END_REF]. WEKA is a comprehensive framework that enables researchers to compare and verify machine learning algorithms.

The organization of the MODBUS data set is similar to that of the KDD Cup 1999 Data Set [START_REF] Hettich | The UCI KDD Archive[END_REF]. Each instance in the data set represents one captured network transaction pair (e.g., merged MODBUS query and response). An instance includes network traffic information and the current state of the process Table 2. Data sets.

Data Set Index

Data Set I Gas pipeline system complete data set Data Set II Water storage tank system complete data set Data Set III Gas pipeline system reduced (10%) data set Data Set IV Water storage tank system reduced (10%) data set control system based on payload content. Note that each instance contains a label identifying it as normal MODBUS traffic or as attack traffic with the designated attack class.

Four data sets were created as part of this research. Table 2 provides the descriptions of the four data sets. Data Set I contains transactions from the gas pipeline system. Data Set II contains transactions from the water storage tank system. The two data sets were generated from network flow records captured with a serial port data logger.

Two reduced size data sets were also created. Data Set III is a gas pipeline system data set, which was created by randomly selecting 10% of the instances in Data Set I. Likewise, Data Set IV is a water storage tank system data set, which was created by randomly selecting 10% of the instances in Data Set II. The two reduced data sets minimize memory requirements and processing time when validating classification algorithms. They are intended for applications for which quick feedback is desired.

Two categories of features are present in the data sets: network traffic features and payload content features. Network traffic features describe the communications patterns in SCADA systems. Compared with traditional enterprise networks, SCADA network topologies and services are relatively static. Note that some attacks against SCADA systems may change network communications patterns. As such, network traffic features are used to describe normal traffic patterns in order to detect malicious activity. Network traffic features include the device address, function code, length of packet, packet error checking information and time intervals between packets. Payload content features describe the current state of the SCADA system; they are useful for detecting attacks that cause devices (e.g., PLCs) to behave abnormally. Payload content features include sensor measurements, supervisory control inputs and distributed control states.

Network Traffic Features

Table 3 lists the ten attributes that comprise the network traffic features. The first and second attributes are the command device address and response device address. Note that the MODBUS serial command address is one byte long, with each server having a unique device address. As such, the command and response device addresses should match during normal operations. An address mismatch is an indicator of a reconnaissance attack. MODBUS serial systems are configured so that all the slave devices (servers) see all the master transactions. Each slave must check the device address to discern the intended recipient before acting on a packet. Based on the system configuration, the set of device addresses that a slave device should encounter is fixed; device addresses not specified in the configuration are anomalous. The command memory, response memory, command memory count and response memory count include internal memory addresses and field sizes for read and write commands. The memory of a MODBUS server is grouped into data blocks called coils, discrete inputs, holding registers and input registers. Coils and discrete inputs represent a single, read-only Boolean bit with authorized values of 0x00 and 0xFF. Holding and input registers are 16-bit words; holding registers are read/write capable while input registers are read only. Each data block may have its own set of contiguous address space or the data blocks may share a common memory address space based on vendor implementation. The command memory and response memory features are coil or register read/write start addresses taken from command and response packets, respectively. The command and response memory count features are the numbers of objects to be read and written, respectively.

The command and response packet length features provide the lengths of the MODBUS query and response frames, respectively. The MODBUS protocol data unit (PDU) is limited to 253 bytes with an additional three bytes for device ID and CRC fields, resulting in a 256-byte packet. In the gas pipeline and water storage tank systems, the master repeatedly performs a block write to a fixed memory address followed by a block read from a fixed memory address. The read and write commands have fixed lengths for each system, and the read and write responses have fixed lengths for each system. Note, however, that many of the described attacks have different packet lengths. As such, the packet length feature provides a means to detect many attacks.

The time interval attribute is a measurement of the time between a MOD-BUS query and its response. The MODBUS protocol is a request-response protocol and the time interval varies only slightly during normal operations. The malicious command injection, malicious response injection and DOS attacks often result in significantly different time interval measurements due to the nature of the attacks. The last attribute is the command/response CRC error rate. This attribute measures the rates of CRC errors identified in command and response packets. Because SCADA network traffic patterns are relatively static, the normal command and response CRC error rates are expected to stay somewhat constant. In a normal system, the error rates should be low; however, the rates are expected to increase when a system is subjected to a denial-of-service attack such as the invalid CRC attack.

Payload Content Features

The payload content features differ for the gas pipeline and water storage system data sets due to different control schemes and different measured variables. The attributes common to both systems are listed in Table 4. During normal operations, the response function code matches the command function code if there is no error. If there is an error, the response sub-function code is the command function code value plus 0x80. The measurement attribute provides the current value of the gas pipeline pressure or water tank level. The naive malicious response injection attack and the complex malicious response injection attack influence process measurements by manipulating the expected values. The system control mode is determined based on data in a command packet. The system control mode can place the system in the shutdown, manual or automatic modes; zero represents the shutdown mode, one represents the manual mode and two represents the automatic mode. A malicious state command injection attack can attempt to modify the system operating mode or shut down the system. The gas pipeline system/water storage tank system use a compressor/pump to add air/water, respectively, to maintain the desired setpoint. If the compressor/pump state has a value of one, then the compressor/pump is on; if it is zero, the compressor/pump is off. When a system is in the automatic mode, the PLC logic controls the compressor/pump state. A malicious complex response injection attack may modify this value in order to mask the actual compressor/pump working state. Note that, in the manual mode, the compressor/pump state is controlled by the manual compressor/pump setting value. A malicious state command injection attack may change the compressor/pump mode continually or intermittently. Table 5 shows the eight attributes that are specific to the gas pipeline system. The initial attribute identifies the setpoint for the nominal gas pressure. The second attribute identifies the operating mode of the system. In the automatic mode, the PLC logic attempts to maintain the gas pressure in the pipeline using a PID control scheme by selecting if the compressor or the relief valve is activated. If the control scheme is zero, then the compressor is activated to increase pressure; if the control scheme is one, then the relief valve is activated using a solenoid to decrease the pressure. In the manual mode, the operator controls the pressure by sending commands to start the compressor or open the relief valve. Additionally, there are five attributes related to the PID controller. The gain, reset, dead band, rate and cycle time impact PID controller behavior and should be fixed during system operation. A malicious parameter command injection attack tries to modify these parameters to interrupt normal control operations. Table 6 shows the four attributes that are specific to the water storage tank system: HH, H, L and LL. In the automatic mode, the PLC logic maintains the water level between the L and H setpoints using an on/off controller scheme. When the sensors detect that the water level has reached the L level, the PLC logic turns the water pump on. Alternatively, when the sensors determine that the water level has reached the H level, the PLC logic turns the water pump off. Note that the water storage tank includes a manual drainage valve that allows water to drain out of the tank when the valve is open. If the manual drainage valve is open, the water level in the tank oscillates between the H and L setpoints continuously as the pump cycles on and off to compensate. When the manual drainage valve is closed, the pump stays on until the water level reaches the H setpoint, at which point it turns off and maintains a constant level. Due to a system fault, if the water level rises to the HH setpoint or falls to the LL setpoint, then an alarm is triggered at the human machine interface that monitors the water storage tank. In the manual mode, the pump state is controlled manually by the human machine interface (i.e., an operator can manually activate and deactivate the pump). Table 7 lists the eight possible label values. Recall that each data set instance is labeled as normal or according to its attack class. The labeling scheme was chosen to match the KDD Cup 1999 Data Set [START_REF] Hettich | The UCI KDD Archive[END_REF], which identified attacks by class. Note that specific attacks in each attack class have similar exploit methods and similar impact on the SCADA system.

Discussion

The data sets described in this paper are relevant to other SCADA systemssystems that use protocols other than MODBUS as well as systems other than gas pipelines and water storage tanks. The features in the data sets are divided into two groups in a similar manner as SCADA protocols divide packets into network traffic related fields and content fields. Indeed, other protocols include similar, albeit not identical, network traffic information such as addresses, function codes, payloads and checksums. Additionally, most SCADA protocols tend to adhere to query-response traffic patterns similar to MODBUS. The content features in the data sets include remote commands and system states similar to how other types of systems monitor and update system settings. As such, the data sets provide a framework to measure the accuracy of intrusion detection approaches designed for a variety of SCADA systems.

Conclusions

Researchers have developed numerous intrusion detection approaches for detecting attacks against SCADA systems. To date, researchers have generally engaged unique threat models and the associated network traffic data sets to train and validate their intrusion detection systems. This leads to a situation where researchers cannot independently verify the results of other research efforts, cannot compare the effectiveness of intrusion detection systems against each other and ultimately cannot adequately judge the quality of intrusion detection systems.

The four data sets developed in this research include network traffic, process control and process measurement features from two laboratory-scale SCADA systems. Data Set I contains transactions from a gas pipeline system while Data Set II contains transactions from a water storage tank system. The data sets were generated from network flow records captured with a serial port data logger in a laboratory environment. A set of 28 attacks was used to create the data sets; the attacks were grouped into four categories: reconnaissance, response injection, command injection and denial-of-service attacks. Reduced size data sets corresponding to Data Sets I and II were also created. Data Set III is a gas pipeline system data set containing 10% of the instances in Data Set I while Data Set IV is a water storage tank system data set containing 10% of the instances in Data Set II. The four data sets comprising normal and attack traffic can be used by security researchers to compare different SCADA intrusion detection approaches and implementations.

Figure 1 .

 1 Figure 1. Gas pipeline and water storage tank systems.

Table 1 .

 1 Intrusion detection systems by threat model and network protocol.

Table 3 .

 3 Attacks on MODBUS systems.

	Attribute	Description
	command address	Device ID in command packet
	response address	Device ID in response packet
	command memory	Memory start position in command packet
	response memory	Memory start position in response packet
	command memory count Number of memory bytes for R/W command
	response memory count	Number of memory bytes for R/W response
	command length	Total length of command packet
	response length	Total length of response packet
	time	Time interval between two packets
	crc rate	CRC error rate

Table 4 .

 4 List of common payload attributes.

	Feature Name	Description
	comm fun	Value of command function code
	response fun	Value of response function code
	sub function	Value of sub-function code in the command/response
	measurement	Pipeline pressure or water level
	control mode	Automatic, manual or shutdown
	pump state	Compressor/pump state
	manual pump setting Manual mode compressor/pump setting
	label	Manual classification of the instance

Table 5 .

 5 Unique features of the gas pipeline system data sets.

	Feature Name Description
	set point	Target pressure in the gas pipeline
	control scheme	Control scheme of the gas pipeline
	solenoid state	State of solenoid used to open the gas relief valve
	gain	Gain parameter value of the PID controller
	reset	Reset parameter value of the PID controller
	dead band	Dead band parameter value of the PID controller
	rate	Rate parameter value of the PID controller
	cycletime	Cycle time parameter value of the PID controller

Table 6 .

 6 Unique features of the water storage system data sets.

	Feature Name Description
	HH	Value of HH setpoint
	H	V a l u e o f H s e t p o i n t
	L	V a l u e o f L s e t p o i n t
	LL	Value of LL setpoint

Table 7 .

 7 Instance classification values.

	Label	Label Label
	Name	Value Description
	Normal	0	Instance is not part of an attack
	NMRI	1	Naive malicious response injection attack
	CMRI	2	Complex malicious response injection attack
	MSCI	3	Malicious state command injection attack
	MPCI	4	Malicious parameter command injection attack
	MFCI	5	Malicious function command injection attack
	DoS	6	Denial-of-service attack
	Reconnaissance	7	Reconnaissance attack


Riccardo Santini

Chiara Foglietta

Stefano Panzieri

EVIDENCE THEORY FOR CYBER-PHYSICAL SYSTEMS

Keywords:

Telecommunications networks are exposed to new vulnerabilities and threats due to interdependencies and links between the cyber and physical layers. Within the cyber-physical framework, data fusion methodologies such as evidence theory are useful for analyzing threats and faults. Unfortunately, the simple analysis of threats and faults can lead to contradictory situations that cannot be resolved by classical models.

Classical evidence theory extensions, such as the Dezert-Smarandache framework, are not well suited to large numbers of hypotheses due to their computational overhead. Therefore, a new approach is required to handle the complexity while minimizing the computational overhead. This paper proposes a hybrid knowledge model for evaluating the intersections among hypotheses. A hybrid frame of discernment is presented using a notional smart grid architecture that transforms the basic probability assignment values from the classical framework. Several analyses and simulations are conducted, with the goal of decreasing conflicts between two independent sources. A comparative analysis is performed using different frames of discernment and rules in order to identify the best knowledge model. Additionally, a computational time analysis is conducted.

Introduction

The pervasive growth of network technology has led to the integration of telecommunications technologies and physical processes to create cyber-physical systems. Cardenas, et al. [START_REF] Cardenas | Secure control: Towards survivable cyber-physical systems[END_REF] define a cyber-physical system as integrating computing, communications and storage capabilities with monitoring and/or control of entities in the physical world, which is done in a dependable, safe, secure and efficient manner under real-time constraints. A cyber-physical system is characterized by the tight connection and coordination between cyber and physical resources. Poovendran [START_REF] Poovendran | Cyber-physical systems: Close encounters between two parallel worlds[END_REF] notes that the concept of a cyber-physical system changes the notion of a physical system to include humans, the infrastructure and the software platform in which the overall system is highly networked.

Examples of cyber-physical systems include supervisory control and data acquisition (SCADA) systems that monitor and control electric power grids, oil and gas pipelines, water supply networks and wastewater treatment systems [START_REF] Burmester | Modeling security in cyber-physical systems[END_REF]. Research activities related to these systems usually focus on reliability and resilience. Krishna and Koren [START_REF] Krishna | Adaptive fault-tolerance for cyber-physical systems[END_REF] have proposed an adaptive control methodology for cyber-physical systems to handle failures of cyber and physical components. Cardenas, et al. [START_REF] Cardenas | Secure control: Towards survivable cyber-physical systems[END_REF] have studied integrity, confidentiality and denial-of-service attacks on cyber-physical systems. This paper considers cyber-physical systems in the context of evidence theory, with the goal of properly identifying the causes of faults and threats when a cyber attack compromises power grid operations. Evidence theory has been applied in multi-sensor fusion problems such as diagnosis [START_REF] Basir | Engine fault diagnosis based on multi-sensor information fusion using Dempster-Shafer evidence theory[END_REF]. Siaterlis and Genge [START_REF] Siaterlis | Theory of evidence-based automated decision making in cyber-physical systems[END_REF] have proposed an evidence theory framework for anomaly detection. In contrast, this paper proposes a hybrid knowledge model for evaluating the intersections among hypotheses. The new approach handles complexity while reducing the computational overhead.

Evidence Theory

Evidence theory is a mathematical formalism for handling uncertainty by combining evidence from different sources to converge to an accepted belief [START_REF] Shafer | A Mathematical Theory of Evidence[END_REF]. The basic concept is to reduce uncertainty in order to identify the set that contains the correct answer to a question.

Frame of Discernment

Let Ω = {ω 1 , • • • , ω n } be the frame of discernment -the set of hypotheses that represents a possible value of the variable ω. In classical evidence theory, the hypotheses are assumed to be mutually exclusive [START_REF] Dempster | Upper and lower probabilities induced by a multivalued mapping[END_REF][START_REF] Shafer | A Mathematical Theory of Evidence[END_REF].

Given a frame of discernment Ω, it is possible to define the power set

Γ(Ω) = {γ 1 , • • • , γ 2 |Ω| } with cardinality |Γ(Ω)| = 2 |Ω| .
This set contains all possible subsets of Ω, including the empty set γ 1 = ∅ and the universal set (frame of discernment) γ 2 |Ω| = Ω.

Basic Probability Assignment

Smets and Kennes [START_REF] Smets | The transferable belief model[END_REF] have defined a model for evidence theory called the transferable belief model. The model relies on a basic probability assignment (BPA) function: m : Γ(Ω) → [0, 1]. The BPA function assigns a value between 0 and 1 to each element of the power set subject to the constraint:

γa⊆Γ(Ω) m(γ a ) = 1 with m(∅) = 0. ( 1
)
Each element γ a with m(γ a ) ̸ = 0 is called a focal set. One of the key goals is to quantify the confidence of propositions of the form: "the true value of ω i is in γ a " where γ a ∈ Γ(Ω). For γ a ∈ Γ(Ω), m(γ a ) is the portion of confidence that supports exactly γ a . This means that the true value is in the set γ a ; however, due to the absence of additional information, it is not possible to better support any strict subset of γ a . Note that this does not correspond to a probability function and it does not respect the property of additivity, i.e., m(γ

a ∪ γ b ) ̸ = m(γ a ) + m(γ b ).
Each BPA is an atomic element in the transferable belief model. In fact, each sensor, agent and node must be able to assign BPA values based on subjective assumptions or using algorithms that automatically determine the assignments.

Combination Rules

In the case of independent information sources, a rule that aggregates the data is required. Several combination rules have been proposed in the literature. The most commonly rules are Dempster's rule [START_REF] Dempster | Upper and lower probabilities induced by a multivalued mapping[END_REF] and Smets' rule [START_REF] Smets | The transferable belief model[END_REF]. This paper considers an additional rule, called proportional conflict redistribution no. 6 (PCR-6), in order to obtain sufficient solutions in terms of a qualityconflict ratio.

Dempster's Rule. Dempster's rule of combination [START_REF] Dempster | Upper and lower probabilities induced by a multivalued mapping[END_REF], which was the first to be formalized, is a purely conjunctive operation. This rule strongly emphasizes the agreement between multiple sources and ignores conflicting evidence through a normalization factor:

Dempster{m i , m j }(∅) = 0 (2) Dempster{m i , m j }(γ a ) = γ b ∩γc=γa m i (γ b )m j (γ c ) 1 - γ b ∩γc=∅ m i (γ b )m j (γ c ) ∀γ a ∈ Γ(Ω). (3)
Note that Dempster's rule assigns a null mass to the empty set, which has certain limitations when the conflict value is very high.

Smets' Rule. Smets' rule of combination [START_REF] Smets | The transferable belief model[END_REF] provides the ability to explicitly express contradictions in the transferable belief model by letting m(∅) ̸ = 0. Smet's rule, unlike Dempster's rule, avoids normalization while preserving commutativity and associativity. The rule is formalized as follows:

Smets{m i , m j }(γ a ) = m i (γ a ) ⊗ m j (γ a ) ∀γ a ∈ Γ(Ω) ( 4 )
where

m i (γ a ) ⊗ m j (γ a ) = γ b ∩γc=γa m i (γ b )m j (γ c ) ∀γ a ∈ Γ(Ω). (5)
The inequality m(∅) > 0 can be explained in two ways. The first is the open world assumption of Dempster [START_REF] Dempster | Upper and lower probabilities induced by a multivalued mapping[END_REF], which expresses the idea that the frame of discernment must contain the true value. Necessarily, if the open world assumption is true, then the set of hypotheses must contain all the possibilities. Under this interpretation, if ∅ is the complement of Ω, then mass m(∅) > 0 represents the case where the truth is not contained in Ω.

The second interpretation of m(∅) > 0 is that there is some underlying conflict between sources. Hence, the mass m(∅) represents the degree of conflict. In particular, the mass m(∅) is computed as:

m i (∅) ⊗ m j (∅) = 1 - γ b ∩γc=∅ (m i (γ b ) ⊗ m j (γ c )) . (6)
PCR-6 Rule. The proportional conflict redistribution rule no. 6 (PCR-6) [START_REF]Advances and Applications of DSmT for Information Fusion (Collected Works)[END_REF] is a non-Bayesian rule for combining BPAs. PCR-6 considers two sources of information evaluated as PCR 6 (∅) = 0 and ∀γ a ∈ Γ(Ω) \ ∅ according to the following equation:

PCR 6 {m i , m j }(γ a ) = Smets{m i , m j }(γ a )+ γ b ∈ Γ(Ω) \ γ a , γ a ∩ γ b = ∅ m 2 i (γ a )m j (γ b ) m i (γ a ) + m j (γ b ) + m 2 j (γ a )m i (γ b ) m j (γ a ) + m i (γ b ) . ( 7
)
The conflict is redistributed between the elements of the power set. In the case of high-conflict sources, only the focal sets that generate the conflict are involved in the redistribution (see the normalization factor in Equation ( 7)). Therefore, the solutions obtained after the combination are better in terms of the quality-conflict ratio.

Architecture for Smart Grid Diagnostics

A smart grid is an excellent example of a cyber-physical system -it comprises the physical electrical grid and an integrated telecommunications network that monitors and controls the energy flow. Figure 1 shows a simplified cyberphysical representation of a smart grid. Note that the EMS/DMS control system uses a telecommunications network to send and receive information from substations in the power grid.

Two assumptions are made about the smart grid architecture. The first assumption concerns the information exchanged by the equipment: under normal conditions, the cyber information can be represented by the timing and volume of four packet types (Command, Ack-Receive, Reply and Ack-Response). The second assumption concerns the sensors used for smart grid management: a packet-sniffing sensor is used in the cyber layer to detect the number of packets in the network and a physical layer sensor is used to indicate whether a piece of equipment (e.g., circuit breaker) is working or not.

In order to apply evidence theory to determine the cause of a malfunction, it is necessary to define the appropriate frame of discernment Ω. In the the example under consideration, there are three hypotheses: normal behavior (N ), physical fault (P ) and cyber threat (C). The system has normal behavior when the breaker is working and the network packets conform to the operational timing and volume constraints. A physical fault exists when the sensors detect a breaker fault. A cyber threat exists when there is excess or low packet volume. As shown in Figure 2, in the classical evidence theory framework, the hypotheses are mutually exclusive with empty intersections.

A plausible scenario is simulated using the specified architecture and parameters. The scenario involves an attacker who compromises the operation of a piece of equipment (circuit breaker) via a telecommunication attacks (distributed denial-of-service attack). A simulation, which has a duration of 100 seconds, is divided into four different situations: Situation 1 (0 to 27 seconds): The smart grid behaves normally and no alarms are detected. The breaker is working and the number of network packets in the specified time window is normal. Situation 2 (28 to 35 seconds): The cyber sensor detects an increasing number of packets in the network (due to the attacker's intrusion), but the breaker is still working.

Situation 3 (36 to 95 seconds): The cyber sensor and the physical sensor both detect anomalous behavior. The packet-sniffing sensor detects a high number of packets and the breaker does not respond to commands.

Situation 4 (96 to 100 seconds): The smart grid is back to normal after the cyber-physical attack because the countermeasures were successful. 1 summarizes the simulation events, with a focus on the time and information sources.

The goal is to fuse all the data provided by the sensors during a simulation in order to detect a cyber-physical attack. As such, the relative frame of discernment Ω according to the classical evidence theory is:

Ω = {C, P, N } . ( 8
)
Starting with Ω, the power set is:

Γ(Ω) = {∅, C, P, N, C ∪ P, C ∪ N, P ∪ N, C ∪ P ∪ N } . (9)
Each sensor has to distribute a unitary mass over specific focal sets during a simulation. Using a combination rule, a fusion result can then be obtained.

Specifically, the focal sets for the cyber sensor are {C, N, P ∪ N, Ω}. Note that a cyber security expert could identify a cyber anomaly, but is unlikely to discern a physical anomaly. Similarly, the focal sets for the physical sensor are {P, N, C ∪ N, Ω}.

Santini, et al. [START_REF] Santini | Evidence theory for smart grid diagnostics[END_REF] have used the PCR-6 rule to develop metrics for identifying the effects of cyber attacks that are designed to inflict physical damage. A cyber-physical fault is detected in the presence of mutually exclusive hypotheses by noticing the existence of non-zero similar masses in the cyber cause set and the physical cause set. Such problems are primarily related to the BPA assignments for the sources, which are application dependent. Another problem relates to the interpretation of conflict values that is done in an ad hoc manner. The following exponential function (depending on the number of captured packets) is used as the BPA assignment to set the mass of {C}:

e -(a•p)/x ( 10
)
where a and p are positive tuning parameters and x is the number of packets. Equation ( 10) is used in the same manner to express the mass of {P } after a physical fault, where x is the persistence of the fault. When two information sources that have high conflict exist in the cyber and physical realms, the rough values obtained after fusion using the PCR-6 rule are unsuitable. The solution proposed in [START_REF] Santini | Evidence theory for smart grid diagnostics[END_REF] is to evaluate at each fusion step the conflict value of the mass distribution over Ω using Smet's rule and compare it with the sum of the two masses in {C} and {P }. The cyber-physical alarm triggering equation is given by: max {mPCR-6(γa)} ∀γa ∈ Ω, if mSmets({∅}) ≤ ρ mPCR-6({C}) + mPCR-6({P }) ≥ mSmets({∅}), if mSmets({∅}) ≥ ρ [START_REF]Advances and Applications of DSmT for Information Fusion (Collected Works)[END_REF] where ρ = 0.7 is a pre-defined threshold for an admissible conflict value. Typically, the decision-making rule in evidence theory is set with the highest BPA value after combining the information from all the sources.

In the smart grid case study, Equation ( 11) is not always valid throughout the simulation: during the cyber-physical anomaly, the decision rule yields different sets for the same events (i.e., initially {C} and then {P }).

As shown in Figure 3, the results are quite interesting. During the simulation, m({C}) and m({P }) converge to the same value even if they belong to two exclusive sets as the classical evidence theory assumes.

Using Equation [START_REF]Advances and Applications of DSmT for Information Fusion (Collected Works)[END_REF], it is possible to transmit to the control center the current state of the system, underlying the occurrence of the cyber-physical attack. Upon analyzing the results, it is possible to confirm that an intersection exists among the sets in the frame of discernment.

Smarandache and Dezert [START_REF]Advances and Applications of DSmT for Information Fusion (Collected Works)[END_REF] have proposed an extended version of evidence theory. The extended theory eliminates the constraint on the exclusivity of hypotheses and explicitly considers intersections among the elements of the power set. Although the theory appears to be useful in our case study, the main problem is the intersection operator. In fact, after defining the frame of discernment Ω, it is necessary to define a special power set called the hyper power set D Ω . The cardinality of D Ω due to the intersection operator follows the Dedekind number sequence: 1, 2, 5, 19, 167, 7580, 7828353, 56130437228687557907787... [START_REF]Advances and Applications of DSmT for Information Fusion (Collected Works)[END_REF][START_REF] Wiedemann | A computation of the eighth Dedekin number[END_REF]. Note that only cases up to n < 7 are tractable with current computing technology. This paper resolves the problem by using a hybrid knowledge model based on classical evidence theory and Dezert-Smarandache theory, which is described in the following section.

Exploring the Frame of Discernment

The computational overhead when using the Dezert-Smarandache theory is extremely high. To address this problem, the initial frame of discernment is modified by considering a hybrid knowledge model between classical evidence theory and Dezert-Smarandache theory. In particular, the intersection of {C} and {P } is explicitly evaluated as in the case of Dezert-Smarandache theory, but in the context of classical evidence theory.

The new frame of discernment, which is shown in Figure 4, is given by:

Ω ′ = {C ′ , P ′ , N, C ∩ P } (12)
where

{C ′ } ∈ Ω ′ is equal to {C} \ {P } in the initial frame of discernment Ω, and
{P ′ } ∈ Ω ′ is {P } \ {C} ∈ Ω.
The intersection {C ∩ P } is added to the frame of discernment because most of the conflict is between the sets {C} and {P }. The new power set is given by:

Γ(Ω ′ ) = {∅, C ′ , P ′ , N, C ∩ P, C ′ ∪ P ′ , C ′ ∪ N, C ′ ∪ (C ∩ P ), P ′ ∪ N, P ′ ∪ (C ∩ P ), N ∪ (C ∩ P ), C ′ ∪ P ′ ∪ N, C ′ ∪ P ′ ∪ (C ∩ P ), (13)
C ′ ∪ N ∪ (C ∩ P ), P ′ ∪ N ∪ (C ∩ P ), Ω ′ }
In the new approach, when the intersection C ∩ P is embedded as another hypothesis in Ω ′ , the cardinality of Γ(Ω ′ ) is 16. In contrast, using the Dezert-Smarandache approach and the Dedekind sequence, the cardinality of |Γ(Ω ′ )| is 19. Of course, it is possible to apply the new approach for a number of elements n ≥ 4 to obtain a hybrid power set with cardinality < D Ω . Table 2. BPA assignment for cyber sensor with the new frame (a = 5, p = 2).

Percentage

Number of Packets

m(C ′ ) 55% m(α) 0 .55 • e -(a•p)/x m(C ∩ P ) 45% m(α) 0 .45 • e -(a•p)/x m(N ) 55% (1 -m(α)) 0.55•(1 -e -(a•p)/x m(P ′ ∪ N ∪ (C ∩ P )) 31.5% (1 -m(α)) 0.315•(1 -e -(a•p)/x ) m(Ω ′ ) 13.5% (1 -m(α)) 0.135•(1 -e -(a•p)/x )
Considering the results obtained in the case study above and the results obtained using the approach presented in [START_REF] Santini | Evidence theory for smart grid diagnostics[END_REF], we selected the function defined in Equation [START_REF] Siaterlis | Theory of evidence-based automated decision making in cyber-physical systems[END_REF] for the BPA assignment. The BPA values for the cyber sensor and physical sensor are summarized in Tables 2 and3, respectively. Note that the only difference is related to the BPA assignment of the focal sets: Table 3. BPA assignment for physical sensor with the new frame (a = 5, p = 2). m(C) is divided into the sets {C ′ } and {C ∩ P } belonging to Ω ′ , as reported in Table 2.

Percentage

m(P ) is divided between m({P ′ }) and to m({C ∩ P }) of Ω ′ , as reported in Table 3.

m({P ∪ N }) is now assigned to m({P ′ ∪ N ∪ (C ∩ P )}) and m({C ′ ∪ N }) to m({C ′ ∪ N ∪ (C ∩ P )}
), as reported in Tables 2 and3.

As discussed above, the BPA assignment is still an open question in the context of evidence theory. Indeed, there is no consensus on how to assign the BPA values. Thus, the BPA functions are selected based on the application. Note that the values reported in Tables 2 and3 were obtained after exhaustive tests on the system.

5.

Hybrid Power Set: Simulations and Results

The hybrid power set was tested by fusing the information using the Dempster and PCR-6 rules. Figures 5(a) and 5(b) show comparisons of the evaluations of the conflict between the information sources. Note that the conflict value in Ω ′ is smaller than Ω and is reduced by approximately 11% during the simulation compared with the original case. When Dempster's rule is used, the values are low and demonstrate contradictory behavior. Note that the set P -C is set P ′ in Ω ′ and C -P is C ′ in Ω ′ . As shown in Figure 6, during the cyber-physical anomaly, the values of m(C) and m(P ) are approximately the same (≃ 0.05). Note that m(C ∩ P ) has a higher value (≃ 0.2), but this is not relevant because the conflict value is high.

Figure 7 shows the values of the singletons after fusion using the PCR-6 rule. Note that the set P -C is set P ′ in Ω ′ and C -P is C ′ in Ω ′ . In this case, the dashed line (i.e., m(C ∩ P )) is greater than the others during the cyber-physical anomaly. Upon examining Figure 6, it is seen that the values of m(C ∩ P ) are comparable with m(C) or m(P ) using Ω instead of Ω ′ as the frame of discernment. Therefore, with the hybrid power set, it is possible to manage the intersection between hypotheses to obtain good results. Using the new frame of discernment and the PCR-6 rule, an operator is able to recognize, with the help of the fusion algorithm, a cyber-physical anomaly represented by C ∩ P . With the hybrid frame of discernment, the results can be analyzed using a classical metric (see Equation ( 11)). Note that throughout the simulation there is one element of the power set with the highest value. As such, an operator does not need any other metrics to trigger a particular event (i.e., cyber-physical anomaly).

For the other elements of the power set Γ(Ω ′ ), the sets represented in RAM. The script was executed 100 times. Table 4 reports the means and the variances. The frame of discernment with fewer elements (i.e., Ω) requires less time on the average than Ω ′ , but the time required has greater variance. Note that the performance would improve if a non-interpreted programming language such as Java or C++ were to be used. Nevertheless, the results are encouraging with regard to the application of evidence theory in real-time environments.

Conclusions

The application of evidence theory to diagnose faults in a cyber-physical system is an important topic in critical infrastructure protection. In certain situations, such as when cyber and physical faults are both present, the classical Dempster-Shafer evidence theory is somewhat restrictive. Therefore, it is necessary to redefine the frame of discernment to better represent the knowledge model due to non-empty intersections between hypotheses. The Dezert-Smarandache model explicitly considers the intersection, but it has a high computational overhead due to the cardinality of the hyper power set. The solution, as presented in this paper, is to use a hybrid knowledge model where the intersection is included in the frame of discernment. The results obtained are encouraging. The conflict value is lower and the situation is described by the singleton set {C ∩ P } as having the highest value among the elements of the hybrid power set during a cyber-physical anomaly.

Our research is currently focusing on generalizing evidence theory using different BPA values. An issue requiring further research is defining BPAs for different cyber attacks that seek to inflict physical damage. Another problem is to manage conflicts and understand the source of inconsistent results. Addi-tionally, it is necessary to study of theoretical properties of the hybrid power set.

Figure 1 .

 1 Figure 1. Cyber-physical representation of a smart grid.

Figure 2 .

 2 Figure 2. Representation of the frame of discernment.

Figure 3 .

 3 Figure 3. Results using the PCR-6 rule for singletons.

Figure 4 .

 4 Figure 4. Representation of the new frame of discernment.

  • e -(a•p)/t 0.055 m(C ∩ P ) 45%i m(β) 0 .45 • e -(a•p)/t 0.045 m(N ) 55% (1 -m(β)) 0.55•(1e -(a•p)/t ) 0.495 m(C ′ ∪ N ∪ (C ∩ P )) 31.5% (1 -m(β)) 0.315•(1e -(a•p)/t ) 0.2835 m(Ω ′ ) 13.5% (1 -m(β)) 0.135•(1e -(a•p)/t ) 0.1215m(N ) has the same value because its intersection with the new set is empty and {N } ∩ {C ∩ P } = ∅.

  Conflict and sum of m(C) and m(P ) in Ω.

  Conflict and sum of m(P -C), m(C -P ) and the intersection m(C ∩ P ) in Ω ′ .

Figure 5 .

 5 Figure 5. BPA trends in the power set Ω and hybrid power set Ω ′ .

  Figure 8 are the only ones with non-zero masses. The values m(C ′ ∪ N ∪ (C ∩ P )) (triangle-marked line) and m(P ′ ∪ N ∪ (C ∩ P )) (dotted line) are the same.

Figure 8 .

 8 Figure 8. Results using the PCR-6 rule for the remaining meaningful elements of Ω ′ .

Table 1 .

 1 Events during the simulation and the associated alarms.

	Time (sec) Events	Detecting Sensor
	0 -27	Normal State	-
	28 -35	Cyber Anomaly	Cyber Sensor
	36 -95	Cyber Anomaly + Physical Fault Cyber + Physical Sensors
	96 -100	Normal State	-

Table

Acknowledgement

This research was partially supported by the 7th Framework Programme of the European Union STREP Project under Grant Agreement 285647 (COCK-PITCI -Cybersecurity of SCADA: Risk Prediction, Analysis and Reaction Tools for Critical Infrastructures (www.cockpitci.eu).

Table 4. Computational times for the power sets Γ(Ω) and Γ(Ω ′ ).

Mean Time Variance

Table 4 shows the computational times of the fusion script for the two frames of discernment Ω and Ω ′ . The script, which was written in Matlab [START_REF]MATLAB version 8.0.0[END_REF], was tested on a laptop with a 2.6 GHz quad-core Intel Core i7 processor and 8 GB


Okan Pala

David Wilson

Russell Bent

Steve Linger

James Arnold

Chapter 12 ACCURACY OF SERVICE AREA ESTIMATION METHODS USED FOR CRITICAL INFRASTRUCTURE RECOVERY

Keywords: Service area estimates, recovery, Thiessen polygons, cellular automata

Electric power, water, natural gas and other utilities are served to consumers via functional sources such as electric power substations, pumps and pipes. Understanding the impact of service outages is vital to decision making in response and recovery efforts. Often, data pertaining to the source-sink relationships between service points and consumers is sensitive or proprietary, and is, therefore, unavailable to external entities. As a result, during emergencies, decision makers often rely on estimates of service areas produced by various methods. This paper, which focuses on electric power, assesses the accuracy of four methods for estimating power substation service areas, namely the standard and weighted versions of Thiessen polygon and cellular automata approaches. Substation locations and their power outputs are used as inputs to the service area calculation methods. Reference data is used to evaluate the accuracy in approximating a power distribution network in a mid-sized U.S. city. Service area estimation methods are surveyed and their performance is evaluated empirically. The results indicate that the performance of the approaches depends on the type of analysis employed. When the desired analysis includes aggregate economic or population predictions, the weighted version of the cellular automata approach has the best performance. However, when the desired analysis involves facility-specific predictions, the weighted Thiessen polygon approach tends to perform the best.

Introduction

Electric power, water, natural gas, telecommunications and other utilities are served to consumers using functional sources (facilities) such as power substations, pumps and pipes, switch controls and cell towers. Each of these sources is related to a geographical service area that includes consumers. Data pertaining to the source-sink relationships between service points and consumers is often sensitive or proprietary and is, therefore, unavailable to external entities. During emergencies, decision makers who do not have access to utility information must rely on estimates of service areas derived by various methods. Decision makers have a strong interest in quantifying the accuracy of critical infrastructure service area estimation methods and developing enhanced estimation techniques [START_REF] Fenwick | Electrical substation service-area estimation using cellular automata: An initial report[END_REF][START_REF] Newton | On the methodology of defining substation spheres of influence within an electric vehicle project framework[END_REF][START_REF] Sulewski | A Geographic Modeling Framework for Assessing Critical Infrastructure Vulnerability: Energy Infrastructure Case Study[END_REF].

This paper assesses the accuracy of four methods that are commonly used to estimate infrastructure impact after a disruptive event. The term "impact" refers to the inability of a utility to provide a service, such as power or gas, due to infrastructure damage. The paper focuses on two types of impacts: (i) aggregate impacts, such as economic activity and the population affected by the outage; and (ii) point data impacts, such as whether specific assets are included in an outage. The methods include Voronoi (Thiessen) polygons, Voronoi (Thiessen) polygons with weights, cellular automata and cellular automata with weights. The methods are compared using a reference model of a power distribution network for a mid-sized U.S. city.

Background

Power, gas, water and other infrastructures serve customers in geographical regions called service areas. Although infrastructure operators have detailed information about the source-sink relationships between their assets, this information is neither organized to facilitate large-scale analyses nor is it documented by public regulatory agencies. In addition, the data is often highly sensitive or proprietary.

Determining service areas in the absence of data has long been a problem, but estimating the service areas accurately is very important in disaster recovery situations [START_REF] Castongia | A Demand-Based Resource Allocation Method for Electrical Substation Service Area Delineation[END_REF][START_REF] Fenwick | Electrical substation service-area estimation using cellular automata: An initial report[END_REF][START_REF] Newton | On the methodology of defining substation spheres of influence within an electric vehicle project framework[END_REF][START_REF] Sulewski | A Geographic Modeling Framework for Assessing Critical Infrastructure Vulnerability: Energy Infrastructure Case Study[END_REF]. Typically, the geographic boundary of a service point is required to estimate the source-sink relationships between serving entities (sources) and served entities (sinks). Increasing the accuracy of the estimates could lead to more efficient recovery. Moreover, understanding the comparative merits of different estimation approaches is necessary to enable decision makers to select the right mitigation and remediation strategies in disaster situations. This paper focuses on Voronoi diagram (Thiessen polygon) and cellular automata estimation approaches.

2.1

Voronoi Diagrams (Thiessen Polygons)

Voronoi diagrams are named after the Russian mathematician Georgy Voronoi, who defined them in 1908. Voronoi diagrams are also called Thiessen polygons after Alfred Thiessen, who in 1911, used the approach to estimate the average rainfall of a region from a set of values recorded at individual stations. Three aspects of Voronoi diagrams (Thiessen polygons) have been investigated over the years: (i) modeling natural phenomena; (ii) investigating geometrical, combinatorial and stochastic properties; and (iii) developing computer-based representations [START_REF] Aurenhammer | Voronoi diagrams: A survey of a fundamental geometric data structure[END_REF].

Thiessen polygons have been used in a variety of ways to present and analyze data. The success of the method comes from its ability to uniformly and systematically partition a geographical region. Given points in a Euclidean plane, a Thiessen diagram divides the plane according to a nearest-neighbor rule, where each point is associated with the region of the plane closest to it [START_REF] Aurenhammer | Voronoi diagrams: A survey of a fundamental geometric data structure[END_REF]. To create the boundaries, straight lines are drawn between all the points; from the mid-point of each line, a perpendicular line is drawn at equal Euclidean distances to each joining point. The Thiessen polygons take shape when the perpendicular lines are trimmed at their intersections with other lines (Figure 1(a)). Interested readers are referred to [START_REF] Okabe | Spatial Tessellations: Concepts and Applications of Voronoi Diagrams[END_REF] for additional details about Thiessen polygons and to [START_REF] Akabane | Optimal geographic allocation of power quality control centers by Voronoi diagram[END_REF][START_REF] Held | Creating electrical distribution boundaries using computational geometry[END_REF][START_REF] Newton | On the methodology of defining substation spheres of influence within an electric vehicle project framework[END_REF][START_REF] Okabe | Generalized network Voronoi diagrams: Concepts, computational methods and applications[END_REF] for details about using the approach to generate critical infrastructure service boundaries.

One drawback of the Thiessen polygon approach is that it assumes that each point is homogenous (as shown in Figure 1(a)). This is generally not the case because each source point provides varying degrees of service. For example, electric power substations have different load outputs and natural gas transportation systems have different pressures and output capacities.

Using weights based on source points can enhance Voronoi-based methods such as the Thiessen polygon approach. A weighted approach creates Thiessen polygons by computing the weighted Euclidean distances [START_REF] Dong | Generating and updating multiplicatively weighted Voronoi diagrams for point, line and polygon features in GIS[END_REF][START_REF] Gahegan | Data structures and algorithms to support interactive spatial analysis using dynamic Voronoi diagrams[END_REF]. The approach assigns smaller service areas to critical infrastructure elements with lower out-puts. This approach is potentially more realistic than an approach that uses Thiessen polygons with equal weights. For example, as shown in Figure 1(b), a 2 MW electric power substation serves a smaller area than neighboring power substations with larger power outputs.

Cellular Automata

Cellular automata are discrete computational systems that comprise finite or denumerable sets of homogeneous, simple cells as part of spatially and temporally discrete grid structures [START_REF] Berto | Cellular automata[END_REF]. They are often used to create mathematical models of complex natural systems that contain large numbers of simple and identical components with local interactions [START_REF] Wolfram | Cellular Automata and Complexity: Collected Papers[END_REF].

A cellular automata is formally defined as a system composed of adjacent cells or sites (usually organized as a regular lattice) that evolves in discrete time steps. Each cell represents an internal state from a finite set of states. The states in the automata are updated in parallel according to a local rule that considers the neighborhood of each cell [START_REF] Chopard | Cellular automata modeling of physical systems[END_REF].

The cellular automata approach originated with digital computing in the late 1940s [START_REF] Neumann | The Computer and the Brain[END_REF][START_REF] Neumann | Papers of John von Neumann on Computers and Computing Theory[END_REF][START_REF] Neumann | Theory of Games and Economic Behavior[END_REF]. However, it was first used in geographical science in the 1970s [START_REF] Barto | Cellular Automata as Models of Natural Systems[END_REF][START_REF] Tobler | Cellular geography[END_REF]. The interest in geographical information technologies in the 1990s led to numerous geographical applications [START_REF] Creutz | Self-organized criticality and cellular automata[END_REF][START_REF] Kuo | Distribution system load estimation and service restoration using a fuzzy set approach[END_REF][START_REF] Toole | Automated utility service area assessment under emergency conditions[END_REF][START_REF] Torrens | Geographic automata systems[END_REF][START_REF] Torrens | Modeling gentrification dynamics: A hybrid approach[END_REF][START_REF] Torrens | Cellular automata and urban simulation: Where do we go from here?[END_REF]. In retrospect, the adoption of cellular automata by the geographical science community was natural because both fields intrinsically rely on proximity, adjacency, distance, spatial configuration, spatial composition and diffusion. Cellular automata also share mathematical and algorithmic structures with remote sensing, relational databases and object-oriented programming [START_REF] Torrens | Cellular automata[END_REF].

Although cellular automata have been applied to a variety of fields, cellular automata techniques were not used for service area calculations until the last decade [START_REF] Fenwick | Electrical substation service-area estimation using cellular automata: An initial report[END_REF][START_REF] Linger | Estimating electrical service areas using GIS and cellular automata[END_REF]. Like Thiessen polygon approaches, cellular automata algorithms can be run with equal weights or weights based on the actual substation loads. Tools that use cellular automata approaches to estimate service and outage areas include the Interdependency Environment for Infrastructure Simulation Systems (IEISS) [START_REF] Bush | Interdependent Energy Infrastructure Simulation System -IEISS Version 2.1 Technical Reference Manual[END_REF], TranSims [START_REF] Bush | Interdependent Energy Infrastructure Simulation System -IEISS Version 2.1 Technical Reference Manual[END_REF][START_REF] Fenwick | Electrical substation service-area estimation using cellular automata: An initial report[END_REF][START_REF] Toole | Automated utility service area assessment under emergency conditions[END_REF] and Water Infrastructure Simulation Environment [START_REF] Mcpherson | The Water Infrastructure Simulation Environment (WISE) Project[END_REF][START_REF] Visarraga | Development of a Java based water distribution simulation capability for infrastructure interdependency analyses[END_REF].

Assessment Methodology

Four algorithms are used to estimate service areas for electric power: (i) Thiessen polygons; (ii) Thiessen polygons with weights based on the electric power substation loads; (iii) cellular automata; and (iv) cellular automata with weights based on the electric power substation loads.

An electric power network in a mid-sized U.S. city comprising roughly 150 substations is used in the evaluation. The reference dataset includes the transmission network, substations, power demand and substation service areas. The reference substation service areas, which are polygonal in shape, were drawn up by an electric power system expert. Economic and population information derived from the 2010 LandScan dataset [START_REF] Bright | LandScan 2011 Global Population Database[END_REF] is incorporated, along with the daytime/nighttime population information from [START_REF] Ching | National Urban Database and Access Portal Tool[END_REF][START_REF] Mcpherson | A day-night population exchange model for better exposure and consequence management assessments[END_REF].

The ESRI suite of GIS tools was used to implement the Thiessen polygon approach. The weighted Thiessen polygons were created using the publiclyavailable ArcGIS extension [START_REF] Dong | Generating and updating multiplicatively weighted Voronoi diagrams for point, line and polygon features in GIS[END_REF]. IEISS [START_REF] Bush | Interdependent Energy Infrastructure Simulation System -IEISS Version 2.1 Technical Reference Manual[END_REF] was used to create the cellular automata and weighted cellular automata polygons; this algorithm grows cells in a raster format starting from each source point (i.e., electric power substation) until it runs out of space or electric power resources.

Aggregated Impacts

Aggregated impacts are used in situations where coarse information about service areas is required. Examples include total population, total economic activity and total area. In these situations, an error in the spatial extent of a service area is acceptable as long as the extent of the area produces the correct values. To perform the comparisons, the daytime population of each polygon associated with a substation is computed using each of the four methods. The results obtained for each method are compared with the actual population associated with the substation in the reference model. In the comparisons, the method with the lowest error is considered to exhibit better performance.

The process is repeated for the nighttime population, economic activity indicators and total area. The economic activity indicators include direct, indirect and induced economic impacts, as well as the economic impact on business and employment. Similarly, the approach that yields the lowest error with respect to the reference dataset is considered to have the best performance. Direct economic impact is based on the types of businesses in a service area. Indirect economic impact is based on the suppliers of commodities in a service area. Induced economic impact is based on the reduction in factor income in a service area. The economic impact on business and employment considers the overall effect on known businesses and employment [START_REF]National Infrastructure Simulation and Analysis Center, FastEcon Tool Summary Report: Fiscal Year 2008[END_REF].

Point Data Impacts

The spatial accuracy of a service area is important for certain types of analyses, such as if an infrastructure outage impacts other infrastructure assets. For example, an asset that depends on electric power from a substation may be unable to function if the substation is out of service. The analysis computes the spatial agreement between the reference service areas and the calculated service areas. Spatial accuracy is evaluated using a point accuracy test. The metric uses 10,000 (10K) points randomly located within a study area.

Figures 2 through 5 show the service areas created by the four methods overlaid on the random points. Point analysis assesses the accuracy of matching critical facilities with their corresponding service source points through service areas. This type of analysis is widely used in land cover classification accuracy assessments [START_REF] Congalton | A review of assessing the accuracy of classifications of remotely sensed data[END_REF][START_REF] Congalton | Assessing the Accuracy of Remotely Sensed Data: Principles and Practices[END_REF]. For each randomly-placed point, the service area to which the point belongs in the reference model is determined; the same deter- mination is performed for the Thiessen polygon, weighted Thiessen polygon, cellular automata and weighted cellular automata approaches (Figure 6). This information is used to create error matrices for evaluating the approaches.

The following equation is used to extract the overall accuracy measure from an error matrix [START_REF] Congalton | Assessing the Accuracy of Remotely Sensed Data: Principles and Practices[END_REF]:

Overall Accuracy = k i=1 n ii k
where k is the number of substations, n is the number of sample points and n ii is a cell along the matrix diagonal corresponding to row i and column i.

Figure 7 shows an error matrix. The rows are the values for a particular method (e.g., weighted cellular automata) and the columns are the reference values. Four matrices are created (one for each method) to assess the accuracy of each calculated dataset compared with the reference dataset. For each sample point, the reference and calculated polygons in which the point falls (i and j, respectively) are determined. If both polygons represent the service area for the same substation, then the ID fields match (i.e., i = j), which causes the n ij cell in the matrix to be incremented by one.

Proximity confidence analyses are also performed to evaluate if the proximity to the source point in a polygon affects the accuracy of an estimate. In the analyses, the distance between the source point and the substation (defined as a serving point in the reference dataset) is measured for each polygon. To classify the points uniformly based on their proximity to the serving source point, the distances are normalized based on the size of the service area polygon that overlays the point for each method. This approach enables a decision maker to quantify the quality of the results based on where a point is located within a service area. Facilities located closer to the service source (i.e., electric power substation) have higher confidence values than those that are further away from the service source. This reduction in confidence can, in fact, be quantified.

As an example, consider two hospitals as point data. The first hospital is located 500 yards away from Substation A and the second hospital is two miles away from Substation B. If the service area sizes are the same for both substations, it is reasonable to compare the hospital to substation distances and to calculate the confidence that the hospitals are correctly associated with the substations. However, if the service area of Substation A is much smaller than that of Substation B, then the distances must be normalized.

Normalization and point classification are based on the distance to the source. Let P (s) be the service area polygon for service point s, A be the area of polygon P and r be the radius of a circle with the same area A as the service area polygon P (s). Furthermore, let i be a randomly-placed point in the agreement zone (i.e., region where the reference data polygon and the polygon produced by the service area estimation method overlap), d be the distance between i and service point s, which is normalized and classified as follows:

(i) Point i is classified in Proximity Class #1 (closest 25%) if d < r/4. (ii) Point i is classified in Proximity Class #2 (25-50%) if r/4 < d < r/2. (iii) Point i is classified in Proximity Class #3 (50-75%) if r/2 < d < 3r/4. (iv) Point i is classified in Proximity Class #4 (farthest 25%) if d < 3r/4.
All the points are classified using the normalized distances and the classifications are used to measure the effect of proximity on the accuracy of point data and to quantify the confidence in a method when reference data is unavailable. We also investigate the ability of the methods to accurately estimate the service-sink relationships when neighboring polygons are also considered. This is performed by creating a lookup table for each method. The lookup table lists all the existing service area polygons for all the methods along with the neighboring polygons. The table is used to recalculate the point accuracy values. For each point that is misplaced in the reference dataset, it is determined if the point is correctly associated with a neighboring polygon. Taking into account the neighboring service area polygons makes it possible to test the accuracy of determining source-sink relationships for critical point locations.

Experimental Results

This section presents the results of the performance analysis of the four approaches, namely the standard and weighted versions of the Thiessen polygon and cellular automata approaches. For the weighted methods, peak energy consumption (in MW) is used for the weights.

Figure 8 shows the results obtained for the four approaches. Each sub-figure displays one service area creation method along with the reference dataset. Figure 8(a) compares the Thiessen polygon approach results with the reference set while Figure 8(b) compares the weighted Thiessen polygon approach results with the reference set. Figures 8(c) and 8(d) show the corresponding results for the cellular automata and weighted cellular automata approaches, respectively.

The first set of results pertains to aggregate statistical accuracy. In particular, the area, population and various economic indicators are compared with the results of the reference service areas.

Table 1 shows the mean differences in the daytime and nighttime populations between the calculated and reference service areas. A smaller value is a better result because the population value produced by the method is closer to the population value produced for the reference service area. For the daytime and nighttime populations, the weighted cellular automata approach yields the best results (smallest differences) compared with the reference data. On the other hand, the cellular automata approach yields results with the highest differences. Similar results were obtained for the sum of differences in populations (Table 2). The weighted cellular automata approach yields the best results. The weighted Thiessen polygon approach yields better results than the standard Thiessen polygon and cellular automata approaches. Tables 3 and4 show the means of the differences in the economic impact for various metrics (direct, indirect, induced, employment and business). In all cases, the difference is the lowest for the weighted cellular automata approach, second lowest for the cellular automata approach and third lowest for the weighted Thiessen polygon approach. The only exception is the economic impact on business (Table 4), for which the weighted Thiessen polygon approach and cellular automata approach swap places. The largest mean difference value is produced by the Thiessen polygon approach. Although the mean difference for the weighted Thiessen polygon approach is larger than that for the cellular automata approach (except for the economic impact on business), the differences are not as notable as the differences for the other categories. Tables 5 and6 show the results corresponding to the sums of the differences; the results have the same trends as in the case of the mean differences. The final aggregate statistic is the total surface area of the polygons. The results of the total surface area comparisons indicate that the average reference polygon area is 1,033 acres. As shown in Table 7, the weighted cellular automata and Thiessen polygon approaches yield polygons that are the closest in size (on average) to the reference polygon sizes. The cellular automata approach yields the least accurate approximation for this metric. For the point accuracy analysis, 10,000 points were selected randomly across the study area and an error matrix was created for each method. The matrices were used to calculate the overlay agreement accuracy. Table 8 shows that the weighted Thiessen polygon approach yields the best overall results (68.9%), followed by the weighted cellular automata approach (59.5%), while the cellular automata approach has the lowest accuracy (52.3%).

As shown in Table 9, the results are nuanced. The weighted cellular automata approach has the highest point accuracy (91%) when points in the closest 25% area of each polygon are considered, followed by the weighted Thiessen polygon approach (86%), the cellular automata approach (85%) and the Thiessen polygon approach (81%). Farther away from the source point, a drop in the accuracy of the unweighted approaches (Thiessen polygon and cellular automata) is observed. The accuracies of the weighted approaches decrease considerably, but they are still higher than the accuracies of the unweighted approaches. It is important to note that the accuracies of all the approaches improve dramatically when neighboring polygons are included. Instead of assigning a point to a single polygon, a point is assigned to a single polygon and a neighboring polygon. This relaxes the analysis to indicate that a point is associated with a source facility from a set of source facilities. The corresponding results are shown in Table 10, where the points are correctly assigned to a set of source facilities more than 95% of the time for all four approaches.

Discussion

Critical infrastructures, such as electric power, natural gas, water and telecommunications, provide vital services to society. In the event of an outage, these services must be restored as soon as possible to bring the situation back to normal and reduce the negative impacts of the outage. Several factors make it difficult for decision makers to assess the impacts of an outage. Critical infrastructure networks are inherently complex and the relationships between network elements as well as those between other networks are not well understood. Outage propagation is complicated to trace, especially in the case of an electric power disruption. In addition, information on source-sink relation-ships is not readily available. Therefore, prioritizing restoration and repair for network elements can be an extremely challenging task.

Moreover, critical infrastructure networks are interconnected and it is often the case that networks depend on other networks to function. For example, an electric network provides power to water pumps, which are part of a water network. Likewise, telecommunications towers and hubs also require electricity to function. Therefore, an electric power network outage can cascade within the network as well as to other networks. The accurate determination of service areas is vital to modeling cross-infrastructure effects. Applying four well-known estimation methods, namely standard and weighted Thiessen polygon and cellular automata approaches, to service area determination for electric power networks yields interesting insights. In general, the weighted cellular automata approach is the best performer while the Thiessen polygon approach has the worst performance. However, for points closest to the boundaries of service areas, the weighted Thiessen polygon approach has the best accuracy.

Visual inspection of the weighted cellular automata polygons compared with the reference dataset polygons provides some insights into the point accuracy results. Two situations lead to the lower accuracy of weighted cellular automata polygons in the point accuracy analysis. The first involves weighted cellular automata polygons at the outer edge of the study area and is an artifact of how the cellular automata algorithm is designed. Cellular automata algorithms favor growth in unconstrained regions and, thus, polygons at the edges tend to grow outward rather than inward, leading to unrealistic results. This behavior can be controlled by introducing boundaries that limit cellular automata growth. The second situation occurs for a few cases in the dataset where the ratio of power output for a specific substation to the total service area in the reference dataset is too large (e.g., when some of the power is provided to an industrial complex). Including substations with large outputs and small area coverage in the reference dataset also contributes to errors.

Finally, cellular automata algorithms incorporate several parameters that must be tuned. This study has used "out of the box" parameters for cellular automata to allow for the least-biased comparisons with Thiessen polygon approaches. However, while parameter tuning can dramatically improve the performance of cellular automata approaches, the tuning is highly specific to the application domain.

Conclusions

Sophisticated modeling and simulation tools are vital to enable decision makers to predict, plan for and respond to complex critical infrastructure service outages [START_REF] Tolone | Critical infrastructure integration modeling and simulation[END_REF][START_REF] Wilson | Recommendation-based geovisualization support for reconstitution in critical infrastructure protection[END_REF]. However, modeling and simulation tools cannot function effectively without adequate, good-quality data. Unfortunately, data pertaining to critical infrastructure assets is highly sensitive and is, therefore, difficult to obtain; detailed data about infrastructure dependencies is even more difficult to obtain.

In the absence of data of adequate quantity and quality, the only feasible solution is to rely on estimation methods to predict the impacts of critical infrastructure service outages on populations, regional economies and other critical infrastructure components. The empirical evaluation of service area estimation techniques described in this paper reveals that the weighted cellular automata and weighted Thiessen polygon approaches produce better estimates than their standard (unweighted) counterparts. Also, the results demonstrate that the weighted cellular automata approach has the best aggregate statistical accuracy while the weighted Thiessen polygon approach has the best point accuracy. However, parameter tuning dramatically improves the performance of the cellular automata approach.

Future research will proceed along three directions. First, other critical infrastructures will be investigated to gain an understanding of the aspects that are unique to critical infrastructures and those that are common between critical infrastructures. Second, other comparison metrics will be developed; for example, substation loads (in MW) could be compared with the expected consumptions by populations and businesses in service areas to assess the accuracy of the computed polygons. Third, formal probability-based methods will be investigated to cope with the error and uncertainty that underlie service area algorithms.

Figure 1 .

 1 Figure 1. (a) Thiessen polygons; (b) Thiessen polygons with weights.

Figure 2 .

 2 Figure 2. Point layer (10K) overlaid with Thiessen polygon layer.

Figure 3 .

 3 Figure 3. Point layer (10K) overlaid with weighted Thiessen polygon layer.

Figure 4 .

 4 Figure 4. Point layer (10K) overlaid with cellular automata polygon layer.

Figure 5 .

 5 Figure 5. Point layer (10K) overlaid with weighted cellular automata polygon layer.

Figure 6 .

 6 Figure 6. Point accuracy assessment data preparation flowchart.

Figure 7 .

 7 Figure 7. Error matrix.

Figure 8 .

 8 Figure 8. Service area polygon examples.

Table 1 .

 1 Mean differences in populations.

	Estimation	Daytime	Nighttime
	Approach	Population Population
	Thiessen Polygons	2,640	3,754
	Weighted Thiessen Polygons	2,554	3,505
	Cellular Automata	3,189	4,847
	Weighted Cellular Automata	1,754	541

Table 2 .

 2 Sum of differences in populations.

	Estimation	Daytime	Nighttime
	Approach	Population Population
	Thiessen Polygons	319K	367K
	Weighted Thiessen Polygons	286K	326K
	Cellular Automata	376K	417K
	Weighted Cellular Automata	79K	24K

Table 3 .

 3 Mean differences in economic impact (direct, indirect and induced).

	Estimation	Direct	Indirect	Induced
	Approach	(dollars) (dollars) (dollars)
	Thiessen Polygons	1.07M	1.54M	2.12M
	Weighted Thiessen Polygons	533K	707K	952K
	Cellular Automata	330K	448K	611K
	Weighted Cellular Automata	11K	8K	16K

Table 4 .

 4 Mean differences in economic impact (employment and business).

	Estimation	Employment Business
	Approach	(dollars)	(dollars)
	Thiessen Polygons	6,200	560
	Weighted Thiessen Polygons	2,700	80
	Cellular Automata	1,900	200
	Weighted Cellular Automata	150	10

Table 5 .

 5 Sum of differences in economic impact (direct, indirect and induced).

	Estimation	Direct	Indirect	Induced
	Approach	(dollars) (dollars) (dollars)
	Thiessen Polygons	59M	85M	116M
	Weighted Thiessen Polygons	25M	33M	44M
	Cellular Automata	16M	22M	31M
	Weighted Cellular Automata	500K	360K	714K

Table 6 .

 6 Sum of differences in economic impact (employment and business).

	Estimation	Employment Business
	Approach	(dollars)	(dollars)
	Thiessen Polygons	340K	31K
	Weighted Thiessen Polygons	127K	4K
	Cellular Automata	96K	11K
	Weighted Cellular Automata	6.4K	0.5K

Table 7 .

 7 Average service area polygon size.

	Estimation	Mean	RMS
	Approach	(acres) (acres)
	Reference Data	1,033	2,241
	Thiessen Polygons	1,054	2,145
	Weighted Thiessen Polygons	1,106	2,536
	Cellular Automata	898	1,822
	Weighted Cellular Automata	921	2,299

Table 8 .

 8 Overall accuracy through point analysis.

	Estimation	Acccuracy
	Approach	(%)
	Thiessen Polygons	54.1
	Weighted Thiessen Polygons	68.9
	Cellular Automata	52.3
	Weighted Cellular Automata	59.5

Table 9 .

 9 Proximity confidence analysis accuracy (%).

	Estimation	25%	25%-50% 50%-75% 75%-100%
	Approach	Area	Area	Area	Area
	Thiessen Polygons	81	65	40	31
	Weighted Thiessen Polygons	86	75	58	50
	Cellular Automata	85	68	38	27
	Weighted Cellular Automata	91	76	54	41

Table 10 .

 10 Point accuracy analysis based on polygon neighborhood relaxation.

	Estimation	Acccuracy
	Approach	(%)
	Thiessen Polygons	96.5
	Weighted Thiessen Polygons	97.4
	Cellular Automata	95.2
	Weighted Cellular Automata	97.9