Delay optimal event detection on ad hoc wireless sensor networks

Delay Optimal Event Detection on Ad Hoc Wireless Sensor Networks arXiv:1105.6065v1 [cs.NI] 30 May 2011 K. Premkumar, Venkata K. Prasanthi M., and Anurag Kumar Dept. of Electrical Communication Engineering, Indian Institute of Science, Bangalore, INDIA email: kprem@ece.iisc.ernet.in, prasanthi.m@gmail.com, anurag@ece.iisc.ernet.in We consider a small extent sensor network for event detection, in which nodes take samples periodically and then contend over a random access network to transmit their measurement packets to the fusion center. We consider two procedures at the fusion center to process the measurements. The Bayesian setting is assumed; i.e., the fusion center has a prior distribution on the change time. In the first procedure, the decision algorithm at the fusion center is network–oblivious and makes a decision only when a complete vector of measurements taken at a sampling instant is available. In the second procedure, the decision algorithm at the fusion center is network–aware and processes measurements as they arrive, but in a time causal order. In this case, the decision statistic depends on the network delays as well, whereas in the network–oblivious case, the decision statistic does not depend on the network delays. This yields a Bayesian change detection problem with a tradeoff between the random network delay and the decision delay; a higher sampling rate reduces the decision delay but increases the random access delay. Under periodic sampling, in the network–oblivious case, the structure of the optimal stopping rule is the same as that without the network, and the optimal change detection delay decouples into the network delay and the optimal decision delay without the network. In the network–aware case, the optimal stopping problem is analysed as a partially observable Markov decision process, in which the states of the queues and delays in the network need to be maintained. A sufficient statistic for decision is found to be the network–state and the posterior probability of change having occurred given the measurements received and the state of the network. The optimal regimes are studied using simulation. Categories and Subject Descriptors: C.2.3 [Computer-Communication Networks]: Network Operations—Network monitoring; I.2.8 [Artificial Intelligence]: Problem Solving, Control Methods, and Search—Control theory General Terms: Algorithms, Design, Performance Additional Key Words and Phrases: Optimal change detection over a network, detection delay, cross–layer design of change detection This is an expanded version of a paper that was presented in IEEE SECON 2006. This work was supported in part by grant number 2900 IT from the Indo-French Center for the Promotion of Advanced Research (IFCPAR), and in part by a project from DRDO, Government of India. Permission to make digital/hard copy of all or part of this material without fee for personal or classroom use provided that the copies are not made or distributed for profit or commercial advantage, the ACM copyright/server notice, the title of the publication, and its date appear, and notice is given that copying is by permission of the ACM, Inc. To copy otherwise, to republish, to post on servers, or to redistribute to lists requires prior specific permission and/or a fee. c 20YY ACM 1529-3785/20YY/0700-0001 $5.00 ACM Transactions on Sensor Networks, Vol. V, No. N, Month 20YY, Pages 1–39. 2 · Fusion Centre Fig. 1. An ad hoc wireless sensor network with a fusion center is shown. The small circles are the sensor nodes (“motes”), and the lines between them indicate wireless links obtained after a self-organization procedure. 1. INTRODUCTION A wireless sensor network is formed by tiny, untethered devices (“motes”) that can sense, compute and communicate. Sensor networks have a wide range of applications such as environment monitoring, detecting events, identifying locations of survivors in building and train disasters, and intrusion detection for defense and security applications. For factory and building automation applications, there is increasing interest in replacing wireline sensor networks with wireless sensor networks, due to the potential reduction in costs of engineering, installation, operations, and maintenance [Honeywell Inc] [ISA]. Event detection is an important task in many sensor network applications. In general, an event is associated with a change in the distribution of a related quantity that can be sensed. For example, the event of a gas leakage at any joints in a pipe causes a change in the distribution of pressure at the joint and hence can be detected with the help of pressure sensors. In this paper, we limit our discussion to the centralized fusion model (see Figure 1), in which each mote, in an event detection network, senses and sends some function of its observations (e.g., quantized samples) to the fusion center at a particular rate. The fusion center, by appropriately processing the sequence of values it receives, makes a decision regarding the state of nature, i.e., it decides whether a change has occurred or not. Our problem is that of minimizing the mean detection delay (the delay between the event occurring and the detection decision at the fusion center) with a bound on the probability of false alarm. We consider a small extent network in which all the sensors have the same coverage, i.e., when the change in distribution occurs it is observed by all the sensors and the statistics of the observations are the same at all the sensors. N sensors synchronously sample their environment at a particular rate. Synchronized operation across sensors is practically possible in networks such as 802.11 WLANs and Zigbee networks since the access point and the PAN coordinator, respectively, transmit beacons that provide all nodes with a time reference. Based on the measurement samples, the nodes send certain values (e.g., quantized samples) to the fusion center. Each value is carried by a packet, which is transmitted using a contention–based multiple access mechanism. Thus, our small extent ACM Transactions on Sensor Networks, Vol. V, No. N, Month 20YY. · 3 sampling instants t0 t1 t2 X1 X2 t3 t4 t5 t6 t7 (1) D 2 vectors of values to be sent by the N sensors (2) D 2 (3) D 2 Fig. 2. The sensors take samples periodically at instants t1i, t2 , · · · , and prepare to send to the h (1) (2) (N) fusion center a vector of values Xh = Xh , Xh , · · · , Xh at th . Each is queued as a packet in the queue of the respective node. Due to multiple access delays, the packets arrive with random (1) (2) (3) delays at the fusion center; for example, for X2 , the delays D2 , D2 , D2 , for the packets from sensors 1, 2 and 3, are shown. network problem is a natural extension of the standard change detection problem (see [Veeravalli 2001] and the references therein) to detection over a random access network. The problem of quickest event detection problem in a large extent network (where the region of interest is much larger than the sensing coverage of any sensor) is considered by us in [Premkumar et al. 2009]. Also, a small extent network can be thought of as a cluster in a large extent network and that the decision maker can be thought of as a cluster head. In this setting, due to the multiple access network delays between the sensor nodes and the fusion center, several possibilities arise. In Figure 2 we show that although the sensors take samples synchronously, due to random access delays the various packets sent by the sensors arrive at the fusion center asynchronously. As shown in the figure, the packets generated due to the samples taken at time t2 (1) (2) (3) arrive at the fusion center with a delay of D2 , D2 , D2 , etc. It can even happen that a packet corresponding to the samples taken at time t3 can arrive before one of the packets generated due to the samples taken at time t2 . Figure 3 depicts a general queueing and decision making architecture in the fusion center. All samples are queued in per–node queues in a sequencer. The way the sequencer releases the packets gives rise to the following three cases, the first two of which we study in this paper. (1) The sequencer queues the samples until all the samples of a “batch” (a batch is the set of samples generated at a sampling instant) are accumulated; it then releases the entire batch to the decision device. The batches arrive to the decision maker in a time sequence order. The decision maker processes the batches without knowledge of the state of the network (i.e., reception times at the fusion center, and the numbers of packets in the various queues). We call this, Network Oblivious Decision Making (NODM). In factory and building automation scenarios, there is a major impetus to replace wireline networks between sensor nodes and controllers. In such applications, the first step could be to retain the fusion algorithm in the controller, while replacing the wireline ACM Transactions on Sensor Networks, Vol. V, No. N, Month 20YY. 4 · fusion centre 1 Sensors and Random Access Network 2 Decision Maker N Sequencer Fig. 3. A conceptual block diagram of the wireless sensor network shown in Figure 1. The fusion center has a sequencing buffer which queues out–of–sequence samples and delivers the samples to the decision maker in time–order, as early as possible, batch–wise or sample–wise. network with a wireless ad hoc network. Indeed, we show that this approach is optimal for NODM, provided the sampling rate is appropriately optimized. (2) The sequencer releases samples only in time–sequence order but does not wait for an entire batch to accumulate. The decision maker processes samples as they arrive. We call this, Network Aware Decision Making (NADM). In NADM, whenever the decision maker receives a sample, it has to roll back its decision statistic to the sampling instant, update the decision statistic with the received sample and then update the decision statistic to the current time slot. The decision maker makes a Bayesian update on the decision statistic even if it does not receive a sample in a slot. Thus, NADM requires a modification in the decision making algorithm in the fusion center. (3) The sequencer does not queue any samples. The fusion center acts on the values from the various sampling instants as they arrive, possibly out of order. The formulation of such a problem would be an interesting topic for future research. Our Contributions: We find that, in the existing literature on sequential change detection problems (see discussion on related literature below), it has been assumed that, at a sampling instant, the samples from all the sensors reach the fusion center instantaneously. As explained above, however, in our problem the delay in detection is not only due to the detection procedure requiring a certain amount of samples to make a decision (which we call decision delay), but also due to the random packet delay in the multiple access network (which we call network delay). We work with a formulation that accounts for both these delays, while limiting ourselves to the particular fusion center behaviours explained in cases (1) and (2) above. In Section 2, we discuss the basic change detection problem and setup the model. In Section 3, we formulate the change detection problem over a random access network in a way that naturally includes the network delay. We show that in the case of NODM, the problem objective decouples into a part involving the network delay and a part involving the optimal decision delay, under the condition that the sampling instants are periodic. Then, in Section 4, we consider the special case of a network with a star topology, i.e., all nodes are one hop away from the fusion center and provide a model for contention in the random access network. ACM Transactions on Sensor Networks, Vol. V, No. N, Month 20YY. · 5 In Section 5, we formulate the NADM problem where we process the samples as they arrive at the fusion center, but in a time causal manner. The out–of–sequence packets are queued in a sequencing buffer and are released to the decision maker as early as possible. We show in the NADM case that the change–detection problem can be modeled as a Partially Observable Markov Decision Process (POMDP). We show that a sufficient statistic for the observations include the network–state (which include the queue lengths of the sequencing buffer, network–delays) and the posterior probability of change having occurred given the measurements received and the network states. As usual, the optimal policy can be characterised via a Bellman equation, which can then be used to derive insights into the structure of the policy. We show that the optimal policy is a threshold on the posterior probability of change and that the threshold, in general, depends on the network state. Finally, in Section 6 we compare, numerically, the mean detection delay performance of NODM and a simple heuristic algorithm motivated by NADM processing. We show the tradeoff between the sampling rate r and the mean detection delay. Also, we show the tradeoff between the number of sensors and the mean detection delay. Related Literature: The basic mathematical formulation in this paper is an extension of the classical problem of sequential change detection in a Bayesian framework. The centralized version of this problem was solved by Shiryaev (see [Shiryaev 1978]). The decentralized version of the problem was introduced by Tenny and Sandell [Tenny and Sandell 1981]. In the decentralized setting, Veeravalli [Veeravalli 2001] provided optimal decision rules for the sensors and the fusion center, in the context of conditionally independent sensor observations and a quasiclassical information structure. For a large network setting, Niu and Varshney [Niu and Varshney 2005] studied a simple hypothesis testing problem and proposed a counting rule based on the number of alarms. They showed that, for a sufficiently large number of sensors, the detection performance of the counting rule is close to that of the optimal rule. In a recent article on anomaly detection in wireless sensor networks [Rajasegarar et al. 2008], Rajasegarar et al. have provided a survey of statistical and machine learning based techniques for detecting various types of anomalies such as sensor faults, security attacks, and intrusions. In [Aldosari and Moura 2004] the authors consider the problem of decentralized binary hypothesis testing, where the sensors quantize the observations and the fusion center makes a binary decision between the two hypotheses. Remark: In the existing literature on the topic of optimal sequential event detection in wireless sensor networks, to the best of our knowledge there has been no prior formulation that incorporates multiple access delay between the sensing nodes and the fusion center. Interestingly, in this paper we introduce, what can be called a cross layer formulation involving sequential decision theory and random access network delays. In particular, we encounter the fork–join queueing model (see, for example, [Baccelli and Makowski 1990]) that arises in distributed computing literature. 2. THE BASIC CHANGE DETECTION PROBLEM In this section, we introduce the model for the basic change detection problem. The notation, we follow, is given here. ACM Transactions on Sensor Networks, Vol. V, No. N, Month 20YY. · 6 Sampling instants Beginning of Slot k t b t b+1 ) k−1 k k+1 Slots Sampling interval of 1/r slots Fig. 4. Time evolves over slots. The length of a slot is assumed to be unity. Thus, slot k represents the interval [k, k + 1) and the beginning of slot k represents the time instant k. Samples are taken periodically every 1/r slots, starting from t1 = 1/r. coarse sampling delay t0 t1 U1 T t2 network delay t3 U2 t~ K U3 t4 U4 t5 U5 t6 t7 U~ K change time detection time without n/w delay ~ T detection time with n/w delay ~ U Fig. 5. Change time and the detection instants with and without network delay are shown. The coarse sampling delay is given by tK − T where tK is the first sampling instant after change, and the network delay is given by UK̃ − tK̃ . • Time is slotted and the slots are indexed by k = 0, 1, 2 . . .. We assume that the length of a slot is unity and that slot k refers to the interval [k, k + 1). Thus, the beginning of slot k indicates the time instant k (see Figure 4). • N sensors are synchronously sampling at the rate r samples/slot, i.e., the sensors make an observation every 1/r slots and send their observations to the fusion center. Thus, for example, if r = 0.1, then a sample is taken by a sensor every 10th slot. We assume that 1/r is an integer. The sampling instants are denoted t1 , t2 , . . . (see Figure 5). Define t0 = 0; note that the first sample is taken at t1 = 1/r. • The vector of network delays of the batch b is denoted by h i (1) (2) (N ) Db = Db , Db , · · · , Db (i) where Db ∈ {1, 2, 3, · · · } is the network delay in slots, of the ith component of (i) the bth batch (sampled at tb = b/r). Also, note that Db > 1, as it requires one time slot for the transmission of a packet to the fusion center after a successful contention. • The state of nature Θ ∈ {0, 1}. 0 represents the state “before change” and ACM Transactions on Sensor Networks, Vol. V, No. N, Month 20YY. · 7 1 represents the state “after change”. It is assumed that the change time T (measured in slots), is geometrically distributed i.e., and, for k > 1, P (T = 0) = ρ P (T = k | T > 0) = p(1 − p)(k−1) . (1) The value of 0 for T accounts for the possibility that the change took place before the observations were made. • The vector of outputs from the sensor devices at the bth batch is denoted by h i (1) (2) (N ) Xb = Xb , Xb , · · · , Xb (i) where Xb ∈ X is the bth output at the ith sensor. Given the state of nature, (i) Xb are assumed to be (conditionally) independent across sensors and i.i.d. over sampling instants with probability distributions F0 (x) and F1 (x) before and after the change respectively. X1 corresponds to the first sample taken. In this work, we do not consider the problem of optimal processing of the sensor measurements to yield the sensor outputs, e.g., optimal quantizers (see [Veeravalli 2001]). • Let Sb , b > 1, be the state of nature at the bth sampling instant and S0 the state at time 0. Then Sb ∈ {0, 1}, with P (S0 = 1) = ρ = 1 − P (S0 = 0) Sb evolves as follows. If Sb = 0 for b > 0, then  1 w.p. pr Sb+1 = 0 w.p. (1 − pr ) where pr = 1 − (1 − p)1/r . Further, if Sb = 1, then Sb+1 = 1. Thus, if S0 = 0, then there is a change from 0 to 1 at the Kth sampling instant, where K is geometrically distributed. For b > 1, P (K = b) = pr (1 − pr )b−1 Each value to be sent to the fusion center by a node is inserted into a packet which is queued for transmission. It is then transmitted to the fusion center by using a contention based multiple access protocol. A node can directly transmit its observation to the fusion center or route it through other nodes in the system. Each packet takes a time slot to transmit. The MAC protocol and the queues evolve over the same time slots. The fusion center makes a decision about the change depending on whether Network Oblivious (NODM) processing or Network Aware (NADM) processing is employed at the fusion center. In the case of NODM processing, the decision sequence (also called as action sequence), is Au , u > 0, with Au ∈ {stop and declare change(1), take another sample(0)}, where u is a time instant at which a complete batch of N samples corresponding to a sampling instant is received by the fusion center. In the case of NADM processing, the decision sequence is Ak , k > 0, with Ak ∈ {stop and declare change(1), take another sample(0)}, i.e., a decision about the change is taken at the beginning of every slot. ACM Transactions on Sensor Networks, Vol. V, No. N, Month 20YY. 8 · fusion centre (2) Xk (N) Xk 1 2 MAC 2 .... X (1) k .... sampling 1 svc. rate N σ Decision Maker N Sequencer Transmitter buffer Fig. 6. A sensor network model of Figure 3 with one hop communication between the sensor nodes and the fusion center. The random access network along with the sequencer is a fork–join queueing model. 3. NETWORK OBLIVIOUS DECISION MAKING (NODM) From Figure 2, we note that although all the components of a batch b are generated (i) at tb = b/r, they reach the fusion center at times tb + Db , i = 1, 2, · · · , N . In NODM processing, the samples, which are successfully transmitted, are queued in a sequencing buffer as they arrive (see Figure 6) and the sequencer releases a (complete) batch to the decision maker, as soon as all the components of a batch arrive. The decision maker makes a decision about the change at the time instants when a (complete) batch arrives at the fusion center. In the Network Oblivious (NODM) processing, the decision maker is oblivious to the network and processes the batch as though it has just been generated (i.e., as if there is no network, hence the name Network Oblivious Decision Making). We further define (see Figure 5) • Ub , (b > 1): the random instant at which the fusion center receives the complete batch Xb e ∈ {0, 1, . . . }: the batch index at which the decision takes place, if there was • K e = 0 means that the decision 1 (stop and declare change) is no network delay. K taken before any batch is generated • Te = tKe : the random time (a sampling instant) at which the fusion center stops and declares change, if there was no network delay e = U e : the random time slot at which the fusion center stops and declares • U K change, in the presence of network delay • Db = Ub − tb : Sojourn time of the bth batch, i.e., the time taken for all the samples of the bth batch to reach the fusion center. Note that Db is given by (i) e at which the max{Db : i = 1, 2, · · · , N }. Thus, the delay of the batch K e − Te detector declares a change is U e − t e = U K K We define the following detection metrics. Mean Detection Delay defined as the expected number of slots between the e in the presence of coarse sampling change point T and the stopping time instant U  h i e −T 1 e and network delays, i.e., Mean Detection Delay = E U {T >T } . ACM Transactions on Sensor Networks, Vol. V, No. N, Month 20YY. · 9 Mean Decision Delay defined as the expected number of slots between the change point T and the stopping time instant Te in the (presence of coarse h samplingdelay i and in the) absence of network delay, i.e., Mean Decision Delay = E Te − T 1{Te>T } . With the above model and assumptions, we pose the following NODM problem: Minimize the mean detection delay with a bound on the probability of false alarm, the decision epochs being the time instants when a complete batch of N components corresponding to a sampling instant is received by the fusion center. In Section 5, we pose the problem of making a decision at every slot based on the samples as they arrive at the fusion center. Motivated by the approach in [Veeravalli 2001] we use the following formulation for a given sampling rate r i h e − T )1 e (2) min E (U {T >T }   such that P Te < T 6 α where α is the constraint on the false alarm probability. ~ T Fig. 7. T ~ U e >T Illustration of an event of false alarm with Te < T , but U Remark 3.1 Note that if α > 1 − ρ, then the decision making procedure can be stopped and an alarm can be raised even before the first observation. Thus, we assume that α < 1 − ρ.   Remark 3.2 Note that here we consider P Te < T as the probability of false   e < T , i.e., a case as shown in Figure 7 is considered a false alarm and not P U alarm. This can be understood as follows: when the decision unit detects a change e , the measurements that triggered this inference would be carrying the at slot U “time stamp” Te, and we infer that the change actually occurred at or before Te. Thus if Te < T , it is an error. We write the problem defined in Eqn. 2 as i h e − T )1 e min E (U {T >T } Πα (3)   where Πα is the set of detection policies for which P Te < T 6 α. Theorem 1 If the sampling is periodic at rate r and the batch sojourn time process Db , b > 1, is stationary with mean d(r), then h i+ i h 1 e −K e − T )1 e min E K = (d(r) + l(r))(1 − α) − ρ · l(r) + min E (U {T >T } Πα r Πα where l(r) is the delay due to (coarse) sampling. ACM Transactions on Sensor Networks, Vol. V, No. N, Month 20YY. 10 · Remark 3.3 For example in Figure 5, the delay due to coarse sampling is t2 − T , e − K = 3 − 2 = 1, and the network delay is U3 − t3 . The stationarity assumption K on Db , b > 1, is justifiable in a network in which measurements are continuously made, but the detection process is started only at certain times, as needed. Proof: The following is a sketch of the proof (the details are in the Appendix – I)  h i h i+  i h e e e e min E (U − T )1{Te>T } = min E (U − T )1{Te>T } + E T − T Πα Πα     h i+  e e + E T −T = min E[D] 1 − P T < T Πα where we have used the fact that under periodic sampling, the queueing system evolution and the evolution of the statistical decision problem are independent, e is independent of {D1 , D2 , . . .} and E[D] is the mean stationary queueing i.e., K delay (of a batch). By writing E[D] = d(r) and using the fact that the problem   i+ h is solved by a policy πα∗ ∈ Πα with P Te < T = α, the problem minΠα E Te − T becomes h i+ h i+ 1 e −K d(r)(1 − α) + min E Te − T = (d(r) + l(r))(1 − α) − ρ · l(r) + min E K Πα r Πα h i+ e − K is the basic where l(r) is the delay due to sampling. Notice that minΠα E K change detection problem at the sampling instants. e and {D1 , D2 , . . .} Remark 3.4 It is important to note that the independence between K arises from periodic sampling. Actually this is conditional independence given the rate of the periodic sampling process. If, in general, one considers a model in which the sampling is at random times (e.g., the sampling process randomly alternates between periodic sampling at two different rates or if adaptive sampling is used) then we can view it as a time varying sampling rate and the asserted independence will not hold. We conclude that the problem defined in Eqn. 2 effectively decouples into the sum h i+ e −K of the optimal delay in the original optimal detection problem, i.e., 1 minΠα E K r as in [Veeravalli 2001], a part that captures the network delay, i.e., d(r)(1 − α), and a part that captures the sampling delay, i.e., l(r)(1 − α) − ρl(r). 4. NETWORK DELAY MODEL From Theorem 1, it is clear that in NODM processing, the optimal decision device and the queueing system are decoupled. Thus, one can employ an optimal sequential change detection procedure (see [Shiryaev 1978]) for any random access network (in between the sensor nodes and the fusion center). Also, NODM is oblivious to the random access network (in between the sensor nodes and the fusion center) and processes a received batch as though it has just been generated. In the case of NADM (which we describe in Section 5), the decision maker processes samples, ACM Transactions on Sensor Networks, Vol. V, No. N, Month 20YY. · 11 Fusion Centre Fig. 8. A sensor network with a star topology with the fusion center at the hub. The sensor nodes use a random access MAC to send their packets to the fusion center. keeping network–delays into account, thus requiring the knowledge of the network dynamics. In this section, we provide a simple model for the random access network, that facilitates the analysis and optimisation of NADM. N sensors form a star topology1 (see Figure 8) ad hoc wireless sensor network with the fusion center as the hub. They synchronously sample their environment at the rate of r samples per slot periodically. At sampling instant tb = b/r, sensor (i) node i generates a packet containing the sample value Xb (or some quantized version of it). This packet is then queued first-in-first-out in the buffer behind the radio link. It is as if each sample is a fork operation that puts a packet into each sensor queue (see Figure 6). The sensor nodes contend for access to the radio channel, and transmit packets when they succeed. The service is modeled as follows. As long as there are packets in any of the queues, successes are assumed to occur at the constant rate of σ (0 < σ < 1) per slot, with the intervals between the successes being i.i.d., geometrically distributed random variables, with mean 1/σ. If, at the time a success occurs, there are n nodes contending (i.e., n queues are nonempty) then the success is ascribed to any one of the n nodes with equal probability. The N packets corresponding to a sample arrive at random times at the fusion center. If the fusion center needs to accumulate all the N components of each sample then it must wait for that component of every sample that is the last to depart its mote. This is a join operation (see Figure 6). It is easily recognized that our service model, in the case of NODM is the discrete time equivalent of generalized processor sharing (GPS – see, for example, [Kumar et al. 2004]), which can be called the FJQ-GPS (fork-join queue (see [Baccelli and Makowski 1990]) with GPS service). In the case of NADM, the service model is just the GPS. In IEEE 802.11 networks and IEEE 802.15.4 networks, if appropriate parameters are used, then the adaptive backoff mechanism can achieve a throughput that is roughly constant over a wide range of n, the number of contending nodes. This 1 Note that Theorem 1 is more general and does not assume a star topology. ACM Transactions on Sensor Networks, Vol. V, No. N, Month 20YY. 12 · η vs n for different transmission rates (C) in RTS−CTS Aggregate throughput η (Mbps) −−−> 7 Analysis Simulation 6 C= 11 Mbps 5 4 C= 5.5 Mbps 3 2 C= 2.2 Mbps 1 0 0 5 10 15 Number of nodes (n)−−−> 20 300 300 250 250 200 200 150 150 η η Fig. 9. The aggregate saturation throughput η of an IEEE 802.11 network plotted against the number of nodes in the network, for various physical layer bit rates: 2.2 Mbps, 5.5 Mbps, and 11 Mbps . The two curves in each plot correspond to an analysis and an NS–2 simulation. 100 100 50 50 0 0 20 40 number of nodes 60 0 0 20 40 number of nodes 60 Fig. 10. The aggregate saturation throughput η of an IEEE 802.15.4 star topology network plotted against the number of nodes in the network. Throughput obtained with default backoff parameters is shown on the left and that obtained with backoff multiplier = 3, is shown on the right. The two curves in each plot correspond to an analysis and an NS–2 simulation. is well known for the CSMA/CA implementation in IEEE 802.11 wireless LANs; see, for example, Figure 9 [Kumar et al. 2008]. For each physical layer rate, the network service rate remains fairly constant with increasing number of nodes. From Figure 10 (taken from [Singh et al. 2008]) it can be seen that with the default backoff parameters, the saturation throughput of a star topology IEEE 802.15.4 network decreases with the number of nodes N , but with the backoff multiplier = 3, the throughput remains almost constant from N = 10 to N = 50 [Singh et al. 2008]; thus, in the latter case our GPS model can be applicable. Theorem 2 The stationary delay D is a proper random variable with finite mean if and only if N r < σ. Proof: See Appendix – II. Thus, for the FJQ–GPS queueing system to be stable, the sampling rate r is chosen σ such that r < N . ACM Transactions on Sensor Networks, Vol. V, No. N, Month 20YY. · 13 5. NETWORK AWARE DECISION MAKING (NADM) In Section 3, we formulated the problem of NODM quickest change detection over a random access network, and showed that (when the decision instants are Uk , as shown in Figure 5) the optimal decision maker is independent of the random access network, under periodic sampling. Hence, the Shiryaev procedure, which is shown to be delay optimal in the classical change–detection problem (see [Shiryaev 1978]), can be employed in the decision device independently of the random access network. It is to be noted that the decision maker in the NODM case, waits for a complete batch of N samples to arrive, to make a decision about the change. Thus, the mean detection delay of the NODM has a network–delay component corresponding to a batch of N samples. In this section, we provide an alternative mechanism of fusion at the decision device called Network Aware Decision Making (NADM), in which the fusion algorithm does not wait for an entire batch to arrive, and processes the samples as soon as they arrive, but in a time–causal manner. We now describe the processing in NADM. Whenever a node (successfully) transmits a sample across the random access network, it is delivered to the decision maker if the decision maker has received all the samples from all the batches generated earlier. Otherwise, the sample is an out–of–sequence sample, and is queued in the sequencing buffer. It follows that, whenever the (successfully) transmitted sample is the last component of the batch that the decision maker is looking for, then the head of line (HOL) components, if any, in the queues of the sequencing buffer are also delivered to the decision maker. This is because, these HOL samples belong to the next batch that the decision maker should process. The decision maker makes a decision about the change at the beginning of every time slot, irrespective of whether it receives a sample or not. In NADM, whenever the decision maker receives a sample, it takes into account the network–delay of the sample while computing the decision statistic. The network–delay is a part of the state of the queueing system which is available to the decision maker. Thus, unlike NODM, the state of the queueing system also plays a role in decision making. In Section 5.1, we define the state of the queueing system. In Section 5.2, we define the dynamical system whose change of state (from 0 to 1) is the subject of interest to us. We define the state of the dynamical system as a tuple that contains the queueing state, the state of nature, and a delayed state of nature. The delayed state of nature is included in the state of the system so that the (delayed) sensor– observations that the decision maker receives at time instant k + 1 depend only on the state, the control, and the noise of the system at time instant k, a property which is desirable to define a sufficient statistic (see page 244, [Bertsekas 2000a]). We explain the evolution of the state of the dynamical system in Section 5.3. In Section 5.4, we formulate the NADM change detection problem and we find a sufficient statistic for the observations in Section 5.5. In Section 5.6, we provide the optimal decision rule for the NADM change detection problem. 5.1 Notation and State of the Queueing System Recall the notation introduced in Section 2. Time progresses in slots, indexed by k = 0, 1, 2 · · · ; the beginning of slot k is the time instant k. Also, the time instant ACM Transactions on Sensor Networks, Vol. V, No. N, Month 20YY. 14 · ∆k t Bk λk tb t b+1 k Slots k+1 Fig. 11. At time k, the decision maker expects samples (or processes samples) from batch Bk . Also, at time k, λk is the number of slots to go for the next sampling instant and ∆k is the number of slots back at which batch Bk is generated. 1 r tB Last component of batch Bk−1 reached at time k tB k−1 k k−1 k ∆ k =0 k+1 ∆ ∆t = 0 =0 k+1 Slots Bk Fig. 12. Illustration of a scenario in which ∆k = 0. If the last component from batch Bk−1 is received at k, and if there is no sampling instant between tBk−1 and k, then ∆k = 0. Also, note in this case that ∆k = ∆k+1 = · · · = ∆tB = 0. In this scenario, at time instants k, k + 1, · · · , tBk , k all the queues at the sensor nodes and at the sequencer are empty, and at time instant tBk +, all sensor node queues have one packet which is generated at tBk . just after the beginning of time slot is denoted by k+ 2 . Recall that the nodes take samples at the instants 1/r, 2/r, 3/r, · · · . We define the state of the queueing system here. Note that the queueing system evolves over slots. • λk ∈ {1, 2, · · · , 1/r} denotes the number of time slots to go for the next sampling instant, at the beginning of time slot k (see Figure 11). Thus, λk :=  1 − k r mod 1 r  . (4) Thus, λ0 = r1 , λ1 = 1r − 1, · · · , and at the sampling instants tb , λtb = 1r . • Bk ∈ {1, 2, 3, · · · } denotes the index of the batch that is expected to be or is being processed by the decision maker at the beginning of time slot k. Note B0 = B1 = · · · = B1/r = 1. Also, note that the batch Bk is generated at time instant Bk /r. • ∆k ∈ {0, 1, 2, · · · } denotes the delay in number of time slots between the time instants k and Bk /r (see Figure 11).   Bk ∆k := max k − ,0 . (5) r 2 Note that the notation t+ denotes a time embedded to the right of t and is different from the notation (x)+ . Recall that (x)+ := max{x, 0}. ACM Transactions on Sensor Networks, Vol. V, No. N, Month 20YY. · A new sample is taken here if (i) (i) λk+1 = 1/r (i) Lk+1 = Lk − 1{Mk =i} + 1{λk+1=1/r} Lk k 15 k+1 slots Departure w.p. Nσ k if (i) Lk > 0 (i) Fig. 13. The evolution of Lk from time slot k to time slot k + 1. If during time slot k, node i transmits (successfully) a packet to the fusion center (i.e., Mk = i), then that packet is flushed out of its queue at the end of time slot k. Also, a new sample is generated (every 1/r slots) exactly (i) at the beginning of a time slot. Thus, Lk+1 , the queue length of sensor node i just after the (i) (i) beginning of time slot k + 1 (i.e., at (k + 1)+) is given by Lk+1 = Lk − 1{Mk =i} + 1{λk+1 =1/r} . Note that the batches of samples taken after Bk /r and up to (including) k are queued either in the sensor node queues or in the sequencing buffer in the fusion center. If at time k, the fusion center receives a sample which is the last sample from batch Bk−1 , then Bk = Bk−1 + 1. If the sampling instant of the Bk th batch is later than k (i.e., Bk /r > k), then ∆k = 0 (up to time Bk /r at which instant, a new batch is generated). This corresponds to the case, when all the samples generated up to time slot k, have already been processed by the decision maker (see Figure 12). In particular, ∆0 = ∆1 = · · · = ∆ r1 −1 = 0. (i) • Lk ∈ {0, 1, 2, · · · } denotes the queue length of the ith sensor node just after the beginning of time slot k (i.e., at time instant k+). The vector of queue (1) (2) (N ) lengths is Lk = [Lk , Lk , · · · , Lk ]. Let Nk ∈ {0, 1, 2, · · · , N } be the number of non–empty queues at the sensor nodes, just after the beginning of time slot k. Nk := N X i=1 1{L(i) >0} k i.e., the number of sensor nodes that contend for slot k is Nk . Hence, using the (i) network model we have provided in Section 4, the evolution of Lk (see Figure 13) is given by the following: (i) L0 (i) Lk+1 = 0     L(i) + 1{λk+1 =1/r}   k (i) = Lk + 1{λk+1 =1/r}      (i) max{Lk − 1, 0} + 1{λk+1 =1/r} w.p. 1 if Nk = 0, w.p. (1 − σ) if Nk > 0, w.p. σ Nk if Nk > 0. Note that when all the samples generated up to time slot k have already been processed by the decision maker and k is not a sampling instant, i.e., ∆k = 0 and λk 6= 1/r, then Lk = 0 (as there are no outstanding samples in the system). For e.g., L1 = L2 = · · · = L1/r−1 = 0. Also, note that just after sampling instant tb , (i) Ltb > 1. ACM Transactions on Sensor Networks, Vol. V, No. N, Month 20YY. 16 · A new sample is received by the fusion center here if (i) Wk (i) Wk+1 k k+1 Mk = j Mk = j slots σ w.p. N k if (j) Lk > 0 (i) Fig. 14. The evolution of Wk from time slot k to time slot k + 1. If a sample from node i is transmitted (successfully) during time slot k (i.e., Mk = i), then it is received by the fusion center at the end of time slot k (i.e. at (k + 1)−). If this sample is from batch Bk , it is passed on to the (i) (i) decision maker. Otherwise, it is queued in the sequencing buffer, in which case Wk+1 = Wk + 1. On the other hand, if a sample from some other node j is transmitted (successfully) during time slot k (i.e., Mk = j 6= i), and if this sample is the last component to be received from batch Bk by the fusion center, then the HOL packet of the ith sequencing queue, if any, is also delivered to (i) (i) (i) the decision maker. Thus, in this case, Wk+1 = max{Wk − 1, 0}. Note that Wk+1 refers to the queue length corresponding to node i at the sequencer, at the beginning of time slot k + 1. • Mk ∈ {0, 1, 2, · · · , N } denotes the index of the node that successfully transmits in slot k. Mk = 0 means that there is no successful transmission in slot k. Thus, from the network model we have provided in Section 4, we note that  if Nk = 0  0 w.p. 1 0 w.p. (1 − σ) if Nk > 0 Mk = (j)  j w.p. σ if Lk > 0, j = 1, 2, · · · , N Nk (i) • Wk ∈ {0, 1, 2, · · · } denotes the queue length of the ith sequencing buffer at time (1) (2) (N ) k. The vector of queue lengths is given by Wk = [Wk , Wk , · · · , Wk ]. Note that Wk = 0 if ∆k = 0, i.e., the sequencing buffer is empty if there are no outstanding samples in the system. In particular, W0 = W1 = · · · = W 1r = 0. (i) The evolution of Wk is explained in Figure 14. If a sample from node i of a batch (i) (i) later than Bk is successfully transmitted during slot k, then Wk+1 = Wk + 1. If a sample from node j of batch Bk is successfully transmitted and if it is the last sample to be received from batch Bk , then the queue lengths of sequencing (i) (i) buffer are decremented by 1, i.e., Wk+1 = max{Wk − 1, 0}. (i) (i) • Rk ∈ {0, 1} denotes whether the sample XBk has been received and processed (i) (i) by the decision maker at time k. Rk = 0 means that the sample XBk has not (i) (i) yet been received by the decision maker and Rk = 1 means that the sample XBk (i) has been received and processed by the decision maker. The vector of Rk s is (1) (2) (N ) (i) (i) given by Rk = [Rk , Rk , · · · , Rk ]. Note that, if Rk = 0, Wk = 0, i.e., the ith sequencing buffer is empty if the sample expected by the decision maker has not yet been transmitted. Also note that when ∆k = 0, Rk = 0, as the samples from the current batch Bk have yet to be generated or have just been generated. (i) (i) We now relate j kthe queue lengths Lk and Wk . Note that at the beginning of k time slot k, 1/r batches have been generated so far, of which Bk − 1 batches are ACM Transactions on Sensor Networks, Vol. V, No. N, Month 20YY. · 17 completely received by the decision maker. In batch Bk , ith sample is received by (i) (i) the decision maker if Rk = 1. Hence, at time k, Bk − 1 + Rk samples generated by node i have been processed by the decision maker and the remaining samples are in the sensor and sequencing queues. Thus, we have   k (i) (i) (i) − (Bk − 1) − Rk (6) Lk + Wk = 1/r   k − Bk /r + 1/r (i) = − Rk 1/r j k (i) ∆k  if k > Bk /r  1/r + 1 − Rk (i) = 1 − Rk if k = Bk /r   (i) −Rk if k < Bk /r Recalling the definition of ∆k , we write the above j k (i) ∆k  + 1 − Rk  1/r (i) (i) Lk + Wk = 1   0 Eqn. as, if ∆k > 0 if ∆k = 0, λk = 1/r if ∆k = 0, λk 6= 1/r. (7) Note that in the above Eqn. ∆k = 0, λk = 1/r (or equivalently k = Bk /r), corresponds to the case when the samples of batch Bk have just been taken and all the samples from all previous batches have been processed. Thus, in this case (i) (i) Lk = 1 (as Wk = 0). In the case of ∆k = 0, λk 6= 1/r (or equivalently k < Bk /r), all the samples from all previous batches have been processed and a new sample (i) (i) from batch Bk is not taken yet. Thus, in this case Lk = 0 (and Wk = 0). Hence, (i) given Qk = [λk , Bk , ∆k , Wk , Rk ], the queue lengths Lk s can be computed as j k (i) (i) ∆k  + 1 − Rk − Wk if ∆k > 0  1/r (i) Lk = φL(i) (Qk ) := 1 if ∆k = 0, λk = 1/r .(8)   0 if ∆k = 0, λk 6= 1/r Also, Nk = φN (Qk ) := N X 1{φL(i) (Qk )>0} . (9) i=1 Thus, the state of the queueing system at time k, can be expressed as Qk = [λk , Bk , ∆k , Wk , Rk ]. Note that the decision maker can observe the state Qk perfectly. The evolution of the queueing system is explained in the next subsection. 5.2 Evolution of the Queueing System The evolution of the queueing system from time k to time k + 1 depends only on Mk , the success/no–success of contention on the random access channel. Note that the evolution of λk is deterministic and that of ∆k depends on Bk . Hence, to describe the evolution of Qk , it is enough to explain the
evolution of Bk , Wk , n and Rk for various cases of Mk . Let Yk+1 ∈ {∅} ∪ ∪N denote the vector of n=1 X samples received, if any, by the decision maker at the beginning of slot k + 1 (i.e., the decision maker can receive a vector of n samples where n ranges from 0 to N ). ACM Transactions on Sensor Networks, Vol. V, No. N, Month 20YY. 18 · At the beginning of time slot k + 1, the following possibilities arise: • No successful transmission: This corresponds to the case i) when all the queues are empty at the sensor nodes (Nk = 0), or ii) when some queues are non–empty at the sensor nodes (Nk > 0), either no queue attempts, or there is more than one attempt (resulting in a collision). In either case, Mk = 0 and the decision maker does not receive any sample, i.e., Yk+1 = ∅. In this case, it is clear that Bk+1 = Bk , Wk+1 = Wk , and Rk+1 = Rk . • Successful transmission of node j’s sample from a later batch: This corresponds to the case, when the decision maker has already received the jth (j) component of the current batch Bk (i.e., Rk = 1) and that it has not received (i) some sample, say i 6= j, from the batch Bk (i.e., Rk = 0, for some i). The received sample (is an out–of–sequence sample and) is queued in the sequencing (j) (j) buffer (Wk+1 = Wk +1). Thus, in this case, Mk = j and the decision maker does not receive any sample, i.e., Yk+1 = ∅. In this case, it is clear that Bk+1 = Bk , Wk+1 = Wk + ej , and Rk+1 = Rk . • Successful transmission of node j’s current sample which is not the last component of the batch Bk : This corresponds to the case when the decision maker has not received the jth component of the batch Bk before time (j) slot k (Rk = 0), and that it has received all the samples that are generated earlier than that of the successful sample. Also, the fusion center is yet to receive P (i) some other component of batch Bk (i.e., N i=1 Rk < N − 1). Thus, in this case, (j) Mk = j and the decision maker receives the sample Yk+1 = XBk . In this case, it is clear that Bk+1 = Bk , Wk+1 = Wk , and Rk+1 = Rk + ej . • Successful transmission of node j’s current sample which is the last component of the batch Bk : This corresponds to the case when the decision maker has not received the jth component of the batch Bk before time slot k (j) (Rk = 0), and that it has received all the samples that are generated earlier than that of the successful sample. Also, this sample is the last component of P (i) batch Bk , that is received by the fusion center. (i.e., N i=1 Rk = N − 1). In this case (along with the received sample), the queues of the sequencing buffer deliver the head of line (HOL) components (which correspond to the batch index Bk + 1), if any, to the decision maker and the queues are decremented by one (i) (i) (Wk+1 = max{Wk − 1, 0}). Thus, Mk = j and the decision maker receives the h i (i′ ) (i′ ) (i′ ) (j) (i) vector of samples Yk+1 = XBk , XBk1+1 , XBk2+1 , · · · , XBkn−1 where Wk > 0 for +1 (i) i ∈ {i′1 , i′2 , · · · i′n−1 }, and Wk = 0 for i ∈ / {i′1 , i′2 , · · · i′n−1 }. In this case, Bk+1 = Bk + 1, Wk+1 = Wk − ei′1 − ei′2 − · · · − ei′n−1 , and Rk+1 = ei′1 + ei′2 + · · · + ei′n−1 . Thus, the state of the queueing system at time k + 1 can be described by Qk+1 = φQ (Qk , Mk ) := [φλ (Qk , Mk ), φB (Qk , Mk ), φ∆ (Qk , Mk ), φW (Qk , Mk ), φR (Qk , Mk )] . In the next subsection, we provide a model of the dynamical system whose state has the state of nature Θk as one of its constituents. The quickest detection of change of Θk from 0 to 1 (at a random time T ) is the focus of this paper. ACM Transactions on Sensor Networks, Vol. V, No. N, Month 20YY. · 19 5.3 System State Evolution Model Let Θk ∈ {0, 1}, k > 0, be the state of nature at the beginning of time slot k. Recall that T is the change point, i.e., for k < T , Θk = 0 and for k > T , Θk = 1, and that the distribution of T is given in Eqn. 1. The state Θk is observed only through the sensor measurements, but these are delayed. We will formulate the optimal NADM change detection problem as a partially observable Markov decision process (POMDP) with the delayed observations. The approach and the terminology used here is in accordance with the stochastic control framework in [Bertsekas 2000a]. At time k, a sample, if any, that the decision maker receives is generated at time Bk /r < k (i.e., samples arrive at the decision maker with a network–delay of ∆k = k − Brk slots). To make an inference about Θk from the sensor measurements, we must consider the vector of states of nature that corresponds to the time instants k − ∆k , k − ∆k + 1, · · · , k. We define the vector of states at time k, Θk := [Θk−∆k , Θk−∆k +1 , · · · , Θk ]. Note that the length of the vector depends on the network–delay ∆k . When ∆k > 0, Θk = [Θ Bk , Θ Bk +1 , · · · , Θk ], and when r r ∆k = 0, Θk is just [Θk ]. Consider the discrete–time system, which at the beginning of time slot k is described by the state Γk = [Qk , Θk ], where we recall that Qk   = λk , Bk , ∆k , Wk , Rk , Θk = [Θk−∆k , Θk−∆k +1 , · · · , Θk ].   Note that Γ0 = r , 1, 0, 0 , Θ0 . At each time slot k, we have the following set of controls {0, 1} where 0 represents “take another sample”, and 1 represents “stop and declare change”. Thus, at time slot k, when the control chosen is 1, the state Γk+1 is given by a terminal absorbing state t; when the control chosen is 0, the state evolution is given by Γk+1 = [Qk+1 , Θk+1 ], where  1 Qk+1 = φQ (Qk , Mk ),   if ∆k+1 = 0   Θk + 1{T =k+1} ,  Θ , Θ , · · · , Θ , Θ + 1 , k k {T =k+1} Θk+1 = h k−∆k k−∆k +1 i if ∆k+1 = ∆k + 1  1  Θ k−∆k + r1 , Θk−∆k + 1r +1 , · · · , Θk , Θk + 1{T =k+1} , if ∆k+1 = ∆k + 1 − r .
(10) =: φΘ Θk , Qk , Mk , 1{T =k+1} where it is easy to observe that Θk + 1{T =k+1} = Θk+1 . When ∆k+1 = ∆k + 1, the batch Bk is still in service, and hence, in addition to the current state Θk+1 = Θk + 1{T =k+1} , we need to keep the states Θk−∆k , Θk−∆k +1 , · · · , Θk . Also, when ∆k+1 = ∆k + 1 − r1 , then the batch index is incremented, and hence, the vector of states that determines the distribution of the observations sampled h at or after Bk+1 /r and before k + 1 is given by Θk−∆k + r1 , Θk−∆k + 1r +1 , · · · , Θk ,  Θk + 1{T =k+1} . Define Ok := 1{T =k+1} , and define Nk := [Mk , Ok ] be the state–noise during time slot k. The distribution of state–noise Nk given the state of the discrete– ACM Transactions on Sensor Networks, Vol. V, No. N, Month 20YY. 20 · 
time system Γk is given by P Mk = m,Ok = oΓk
= [q, θ] and is the product
of the distribution functions, P Mk = mΓk = [q, θ] and P Ok = oΓk = [q, θ] . These distribution functions are provided in Appendix – III. The problem is to detect the change in the state Θk as early as possible by sequentially observing the samples at the decision maker. 5.4 The NADM Change Detection Problem We now formulate the NADM change–detection problem in which the observations from the sensor nodes are sent over a random access network to the fusion center and the fusion center processes the samples in the NADM mode. In Section 5.3, we defined the state Γk = [Qk , Θk ] on which we formulate the NADM change detection problem as a POMDP. Recall that at the beginning of slot k, the decision maker receives a vector of sensor measurements Yk and observes the state Qk of the queueing system. Thus, at time k, Zk = [Qk , Yk ] is the observation of the decision maker about the state of the dynamical system Γk . Let Ak ∈ {0, 1} be the control (or action) chosen by the decision maker after having observed Zk at k. Recall that 0 represents “take another sample”  and 1 represents the action “stop and declare change”. Let Ik = Z[0:k] , A[0:k−1] be the information vector3 that is available to the decision maker, at the beginning of time slot k. Let τ be a stopping time with respect to the sequence of random variables I1 , I2 , · · · . Note that Ak = 0 for k < τ and Ak = 1 for k > τ . We are interested in obtaining a stopping time τ (with respect to the sequence I1 , I2 , · · · ) that minimizes the mean detection delay subject to a constraint on the probability of false alarm.   (11) min E (τ − T )+ such that P (τ < T ) 6 α Note that in the case of NADM, at any time k, a decision about the change is made based on the information Ik (which includes the batch index we are processing and the delays). Thus, in the case of NADM, false alarm is defined as the event {τ < T } and, hence, τ > T is not classified as a false alarm even if it is due to pre–change measurements only. However, in the case of NODM, this is classified as a false alarm as the decision about the change is based on the batches received until time k. Let c be the cost per unit delay in detection. We are interested in obtaining a stopping time τ ∗ that minimizes the expected cost (Bayesian risk) given by   C(c, τ ∗ ) = min E 1{Θτ =0} + c · (τ − T )+ τ " # τX −1 = min E 1{Θτ =0} + c · 1{Θk =1} τ " k=0 = min E gτ (Γτ , Aτ ) + τ = min E τ 3 The "∞ X k=0 τ −1 X # gk (Γk , Ak ) k=0 # gk (Γk , Ak ) notation Z[k1 :k2 ] := Zk1 , Zk1 +1 , · · · , Zk2 ACM Transactions on Sensor Networks, Vol. V, No. N, Month 20YY. (12) · 21 where, as defined earlier, Γk = [Qk , Θk ]. Let θ = [θδ , θδ−1 , · · · , θ1 , θ0 ]. We define for k 6 τ  0,    1, gk ([q, θ], a) = c,    0, if if if if θ0 θ0 θ0 θ0 = 0, a = 0 = 0, a = 1 = 1, a = 0 = 1, a = 1 (13) and for k > τ , gk (·, ·) := 0. Recall that Ak = 0 for k < τ and Ak = 1 for k > τ . Note that Ak , the control at time slot k, depends only on Ik . Thus, every stopping time τ , corresponds to a policy µ = (µ0 , µ1 , · · · ) such that Ak = µk (Ik ), with Ak = 0 for k < τ and Ak = 1 for k > τ . Thus, Eqn. 12 can be written as ∗ C(c, τ ) = min E µ = min µ = min µ " ∞ X # gk (Γk , Ak ) k=0 ∞ X k=0 ∞ X E[gk (Γk , Ak )] (by monotone convergence theorem) E[gk (Γk , µk (Ik ))] (14) k=0 Since Θk is observed only through Ik , we look at a sufficient statistic for Ik in the next subsection. 5.5 Sufficient Statistic In Section 5.2, we have illustrated the evolution of the queueing system Qk and we have shown in different scenarios, the vector Yk received by the decision maker. Recall from Section 5.2 that  ∅, if Mk = 0,    (j)   ∅, if Mk = j > 0, Rk = 1,    N  P  (i) (j)  if Mk = j > 0, Rk = 0, Rk < N − 1  Yk+1,0 , i=1 Yk+1 = N P  (j)   [Yk+1,0 , Yk+1,1 , · · · , Yk+1,n ] , if Mk = j > 0, Rk = 0, Rk(i) = N − 1,    i=1   N  P   1{W (i) >0} = n.  i=1 k (M ) Note that Yk+1,0 corresponds to XBk k . The last part of the above equation corresponds to the last pending sample of the batch Bk arriving at the decision maker at time k + 1, with some samples from batch Bk + 1 (= Bk+1 ) also being released by the sequencer. In this case, the state of nature at the sampling instant of the batch Bk+1 = Bk + 1 is Θk−∆k +1/r . Note that Θk−∆k +1/r is a component of the vector Θk as k − ∆k + 1/r = (Bk + 1)/r < k. Thus, the distribution of ACM Transactions on Sensor Networks, Vol. V, No. N, Month 20YY. · 22 Yk+1,0 , Yk+1,1 , · · · , Yk+1,n is given by  f0 (·), if Θk−∆k = 0 fYk+1,0 (·) = f1 (·), if Θk−∆k = 1 and  f0 (·), if Θk−∆k +1/r = 0 fYk+1,i (·) = , i = 1, 2, · · · , n. f1 (·), if Θk−∆k +1/r = 1 Thus, at time k + 1, the current observation Yk+1 depends only on the previous state Γk , previous action Ak , and the noise of the system Nk . Thus, a  previous 
 sufficient statistic is P Γk = [q, θ] Ik [q,θ]∈S (see page 244, [Bertsekas 2000a]) where S is the set of all states of the dynamical system defined in Sec. 5.3. Let q = [λ, b, δ, w, r]. Note that 
P Γk = [q, θ] Ik 
= P Γk = [q, θ] Ik−1 , Qk , Yk 
= 1{Qk =q} · P Θk = θ Ik−1 , Qk = q, Yk = 1{Qk =q} 
·P [Θk−δ , Θk−δ+1 , · · · , Θk−1 , Θk ] = [θδ , θδ−1 , · · · , θ1 , θ0 ] Ik−1 , Qk = q, Yk 
= 1{Qk =q} · P Θk−δ = θδ Ik−1 , Qk = q, Yk · δ Y j=1 
P Θk−δ+j = θδ−j Θk−δ+j ′ = θδ−j ′ , j ′ = 0, 1, · · · , j − 1, Ik−1 , Qk = q, Yk (15) Observe that 
P Θk−δ+j = θδ−j Θ[k−δ:k−δ+j−2] , Θk−δ+j−1 = 0, Ik−1 , Qk = q, Yk  1 − p, if θδ−j = 0 = p, if θδ−j = 1 and 
P Θk−δ+j = θδ−j Θ[k−δ:k−δ+j−2] , Θk−δ+j−1 = 1, Ik−1 , Qk = q, Yk  0, if θδ−j = 0 = 1, if θδ−j = 1. This is because given Θk−δ , the events {Θk−δ+j = θδ−j }, {Ik−1 , Qk = q, Yk } are conditionally independent. Thus, Eqn. 15 can be written as 
P Γk = [q, θ] Ik 
 Ik−1 , Qk = q, Yk , if θ = 1 1{Qk =q} · P Θk−δ = 1 
 1 δ−j−1  p, if θ = [0, · · · , 0, |{z} 1 , · · · , 1] {Qk =q} · P Θk−δ = 0 Ik−1 , Qk = q, Yk · (1 − p) =  θ j  
 if θ = 0 1{Qk =q} · P Θk−δ = 0Ik−1 , Qk = q, Yk · (1 − p)δ , (16) ACM Transactions on Sensor Networks, Vol. V, No. N, Month 20YY. · e k := Θk−∆ , and define Define Θ k    e k = 1Ik−1 , Qk = [λ, b, δ, w, r], Yk Ψk := P Θ 
= P Θk−δ = 1Ik−1 , Qk = [λ, b, δ, w, r], Yk 
Πk := P Θk = 1Ik−1 , Qk = [λ, b, δ, w, r], Yk 
= P T 6 k Ik−1 , Qk = [λ, b, δ, w, r], Yk . 23 (17) Thus, Eqn. 16 can be written as 
P Γk = [[λ, b, δ, w, r], θ] Ik  1{Q =[λ,b,δ,w,r]} · Ψk , if θ = 1   1 k δ−j−1 · (1 − Ψ ) · (1 − p) p, if θ = [0, · · · , 0, |{z} 1 , · · · , 1] k {Qk =[λ,b,δ,w,r]} (18) =  θj   1{Qk =[λ,b,δ,w,r]} · (1 − Ψk ) · (1 − p)δ , if θ = 0 We now find a relation between Πk and Ψk in the following Lemma. Lemma 1 The relation between the conditional probabilities Πk and Ψk is given by
(19) Πk = Ψk + (1 − Ψk ) 1 − (1 − p)δ Proof. See Appendix – IV. From Eqn. 18 and Lemma 1, it is clear that a sufficient statistic for Ik is νk = [Qk , Πk ]. Also, we show in Appendix – V that νk can be computed recursively, i.e., when Ak = 0, νk+1 = [Qk+1 , Πk+1 ] = [Qk+1 , φΠ (νk , Zk+1 )], and when Ak = 1, νk+1 = t, a terminal state. Thus, νk can be thought of as entering into a terminating (absorbing) state t at τ (i.e., νk = [Qk , Πk ] for k < τ and νk = t for k > τ ). Since νk is sufficient, for every policy µk there corresponds a policy µ ek such that µk (Ik ) = µ ek (νk ) (see page 244, [Bertsekas 2000a]). 5.6 Optimal Stopping Time τ Let Q be the set of all possible states of the queueing system, Qk . Thus the state space of the sufficient statistic is N = (Q × [0, 1]) ∪ {t}. Recall that the action space is A = {0, 1}. Define the one–stage cost function ge : N × A → R+ as follows. Let ν ∈ N be a state of the system and let a ∈ A be a control. Then,  if ν = t 0 c · π if ν = [q, π], a = 0 ge(ν, a) =  1 − π if ν = [q, π], a = 1. Note from Eqn. 13 for k 6 τ that E[gk (Θk , Ak )] = E[gk (Θk , µk (Ik ))]      = E E gk (Θk , µk (Ik )) Ik = E[e g(νk , µ ek (νk ))] ACM Transactions on Sensor Networks, Vol. V, No. N, Month 20YY. 24 · and for k > τ , E[gk (Θk , Ak )] = 0 = E[e g(t, ·)] Since, {νk } is a controlled Markov process, and the one–stage cost function e g(·,
·), the transition probability kernel for Ak = 1 and for Ak = 0 (i.e., P Zk+1 νk ), do not depend on time k, and the optimization problem defined in Eqn. 14 is over infinite horizon, it is sufficient to look for an optimal policy in the space of stationary Markov policies (see page 83, [Bertsekas 2000b]). Thus, the optimization problem defined in Eqn. 14 can be written as C(c, τ ∗ ) = min µ e ∞ X 
 E ge νk , µ ek (νk ) k=0 ∞ X 
 = E ge νk , µ e∗ (νk ) . (20) k=0 Thus, the optimal total cost is given by    ∞ X
 ∗ ∗  J ([q0 , π0 ]) = E e g νk , µ e (νk ) ν0 = [q0 , π0 ] . (21) k=0 The solution to the above problem is obtained following the Bellman’s equation,      ∗ ∗  J ([q, π]) := min 1 − π, cπ + E J (Qk+1 , φΠ (νk , Zk+1 )) νk = [q, π] . (22) where the function φΠ (νk , Zk+1 ) is provided in Appendix – V. Remark 5.1 The optimal stationary Markov policy (i.e., the optimum stopping rule τ ) in general depends on Q. Hence, the decision delay and the queueing delay are coupled, unlike in the NODM case. We characterize the optimal policy in the following theorem. Theorem 3 The optimal stopping rule τ ∗ is a network–state dependent threshold rule on the a posteriori probability Πk , i.e., there exist thresholds γ(q) such that τ = inf{k > 0 : Πk > γ(Qk )} (23) Proof. See Appendix–VI. In general, the thresholds γ(Qk )s (i.e., optimum policy) can be numerically obtained by solving Eqn. 22 using value iteration method (see pp. 88–90, [Bertsekas 2000b]). However, computing the optimal policy for the NADM procedure is hard as the state space is huge even for moderate values of N . Hence, we resort to a suboptimal policy based on the following threshold rule, which is motivated by the structure of the optimal policy. τ = inf{k > 0 : Πk > γ} ACM Transactions on Sensor Networks, Vol. V, No. N, Month 20YY. (24) · 25 where γ is chosen such that P (τ < T ) = α is met. Thus, we have formulated a sequential change detection problem when the sensor observations are sent to the decision maker over a random access network, and the fusion center processes the samples in the NADM mode. The information for decision making now needs to include the network state Qk (in addition to the samples received by the decision maker); we have shown that [Qk , Πk ] is sufficient for the information history Ik . Also, we have provided the structure for the optimal policy. Since, obtaining the optimal policy is computationally hard, we gave a simple threshold based policy, which is motivated by the structure of the optimal policy. 6. NUMERICAL RESULTS Minimizing the mean detection delay not only requires an optimal decision rule at the fusion center but also involves choosing the optimal values of the sampling rate r, and the number of sensors N . To explore this, we obtain the minimum decision delay for each value of the sampling rate r numerically, and the network delay via simulation. 6.1 Optimal Sampling Rate Consider a sensor network with N nodes. For a given probability of false alarm, the decision delay (detection delay without the network–delay component) decreases with increase in sampling rate. This is due to the increase in the number of samples that the fusion center receives within a given time. But, as the sampling rate increases, the network delay increases due to the increased packet communication load in the network. Therefore it is natural to expect the existence of a sampling rate r∗ , with r∗ < σ/N , (the sampling rate should be less than σ/N , for the queues to be stable; see Theorem 2) that optimizes the tradeoff between these two components of detection delay. Such an r∗ , in the case of NODM can be obtained by minimizing the following expression over r (recall Theorem 1). h i+ 1 e −K (d(r) + l(r)) (1 − α) − ρ · l(r) + min E K r Πα h i+ e −K Note that in the above expression, the delay term minΠα E K also depends on the sampling rate r via the probability of change pr = 1 − (1 − p)(1/r) . The delay due to coarse sampling l(r)(1−α)−ρ·l(r) can be found analytically (see Appendix – h i+ e − K by the asymptotic (as α → 0) I). We can approximate the delay minΠα E K delay as N I(f1 ,f|0ln(α)| )+| ln(1−pr )| where I(f1 , f0 ) is the Kullback–Leibler (KL) divergence between the pdfs f1 and f0 (see [Tartakovsky and Veeravalli 2005]). But, obtaining the network–delay (i.e., d(r)(1 − α)) analytically is hard, and hence an analytical characterisation of r∗ appears intractable. Hence, we have resorted to numerical evaluation. The distribution of sensor observations are taken to be N (0, 1) and N (1, 1)4 , before and after the change respectively for all the 10 nodes. We choose the prob4 As usual, N (a, v) denotes a normal distribution with mean a and variance v ACM Transactions on Sensor Networks, Vol. V, No. N, Month 20YY. 26 · Mean Detection Delay (slots) 220 200 NODM: Simulation NODM: Approximation NADM: Simulation 180 160 140 120 100 80 60 0.01 0.015 0.02 0.025 0.03 0.035 Sampling rate, r (samples per slot) Fig. 15. Mean detection delay for N = 10 nodes is plotted against the sampling rate r for both NODM and NADM (defined in Eqn. 24). For NODM, an approximate analysis is also plotted. This was obtained with the prior probability ρ = 0, p = 0.0005, probability of false alarm target α = 0.01, σ = 0.3636 and with the sensor observations being N (0, 1) and N (1, 1), before and after the change respectively. ability of occurrence of change in a slot to be p = 0.0005, i.e., the mean time h i+ e −K until change is 2000 slots. minΠα E K and d(r) are obtained from simulation for α = 0.01 and σ = 0.3636 and the expression for mean detection delay (displayed above) is plotted against r in Figure 15. Note that both NODM and NADM are threshold based, and we obtain the corresponding thresholds for a target PFA = 0.01 by simulation. These thresholds are then used to obtain the mean detection delay by simulation. In Figure 15, we also plot the approximate mean detection delay which is obtained through the expression for l(r) and the approxh i+ | ln(α)| e −K . We study this approximation imation, minΠα E K ≈ N I(f1 ,f0 )+| ln(1−pr )| as this provides an (approximate) explicit expression for the mean decision delay. The delay in the FJQ–GPS does not have a closed form expression. Hence, we still need simulation for the delay due to queueing network. It is to be noted that at k = 0, the size of all the queues is set to 0. The mean detection delay due to the procedure defined in Eqn. 24 is also plotted in Figure 15. As would have been expected, we see from Figure 15 that the NADM procedure has a better mean detection delay performance than the NODM procedure. Note that σ/N = 0.03636 and hence for the queues to be stable (see Theorem 2), the sampling rate has to be less that σ/N = 0.03636 (1/28 < 0.03636 < 1/27). As the sampling rate r increases to 1/28 (the maximum allowed sampling rate), the queueing delay increases rapidly. This is evident from Figure 15. Also, we see from Figure 15 that operating at a sampling rate around 1/34(≈ 0.0294) samples/slot would be optimum. The optimal sampling rate is found to be approximately the same for NODM and NADM. At the optimal sampling rate the mean detection delay of NODM is 90 slots and that of NADM is 73 slots. ACM Transactions on Sensor Networks, Vol. V, No. N, Month 20YY. Mean Decision Delay (slots) · 62 60 60 50 Delay (slots) 58 56 54 52 27 Coarse sampling delay Decision maker delay 40 30 20 10 50 48 0 0 5 10 15 20 25 30 35 40 0 5 10 Number of Sensors, N 15 20 25 30 35 40 Number of Sensors, N Fig. 16. Mean decision delay of NODM procedure for N × r = 1/3 is plotted against the the number of nodes N . The plot is obtained with ρ = 0, p = 0.0005, α = 0.01 and with the sensor observations being N (0, 1) and N (1, 1), before and after the change respectively. The components of the mean decision delay, i.e., the coarse sampling delay (1 − α)l(r) − ρl(r), and the decision i+ h e −K are shown on the right. maker delay, 1 minΠ E K r α 6.2 Optimal Number of Sensor Nodes (Fixed Observation Rate) Now let us consider fixing N × r. This is the number of observations the fusion center receives per slot in a network with N nodes sampling at a rate r (samples per slot). It is also a measure of the energy spent by the network per slot. Since it has been assumed that the observations are conditionally independent and identically distributed across the sensors and over time, it is natural to ask how beneficial it is to have more nodes sampling at a lower rate, when compared to fewer nodes sampling at a higher rate with the number of observations per slot being the same. With p = 0.0005, α = 0.01, and σ = 0.3636, and f0 ∼ N (0, 1) and f1 ∼ N (1, 1), we present simulation results for two examples, the first one being N r = 1/3 (the case of heavily loaded network) and the second one being N r = 1/100 (the case of lightly loaded network, N r ≪ σ). h i+ e −K Figure 16 shows the plot of mean decision delay, l(r)(1−α−ρ)+ 1 minΠα E K r versus the number of sensors when N r = 1/3. As N increases, the sampling rate r = 1/(3N ) decreases and hence the coarse sampling delay l(r)(1−α) increases; this can be seem to be approximately linear by analysis of the expression for l(r) given in Appendix – I. Also, as N increases, the decision maker gets more samples at the deh i+ e −K cision instants and hence the delay due to the decision maker 1 minΠα E K r decreases (this is evident from the right side of Figure 16). Figure 16 shows that in the region where N is large (i.e., N > 20) or N is very small (i.e., N < 5), as N increases, the mean decision delay increases. This is because in this region as N increases, the decrease in the delay due to decision maker is smaller compared to the increase in the delay due to coarse sampling. However, in the region where N is moderate (i.e., for 5 6 N < 20), as N increases, the decrease in the delay due to decision maker is large compared to the increase in the delay due to coarse sampling. Hence in this region, the mean decision delay decreases with N . Therefore, we infer that when N × r = 31 , deploying 20 nodes sampling at 1/60 samples per slot is optimal, when there is no network delay. ACM Transactions on Sensor Networks, Vol. V, No. N, Month 20YY. 28 · Mean Detection Delay (slots) 220 NODM: Simulation NADM: Simulation 200 180 160 140 120 100 80 60 0 5 10 15 20 25 30 35 40 Number of Sensors, N Fig. 17. Mean detection delay for N × r = 1/3 is plotted against the number of nodes N . This was obtained with ρ = 0, p = 0.0005, α = 0.01 σ = 0.3636 and with the sensor observations being N (0, 1) and N (1, 1), before and after the change respectively. Figure 17 shows the mean detection delay (i.e., the network delay plus the decision delay shown in Figure 16) versus the number of nodes N for a fixed N × r = 1/3. As the the number of nodes N increases, the sampling rate r = 1/(3N ) decreases. For large N (and equivalently small r), in the case of NODM with the Shiryaev procedure, the network delay, d(r) ≈ N σ as it requires N (independent) successes, each with probability σ, in the random access network to transport a batch of N samples (also, since the sampling rate r is small, one would expect that a batch is delivered before a new batch is generated) and the decision maker requires just one batch of N samples to stop (after the change occurs). Hence, for large N , the detection delay is ≈ l(r)(1 − α) + d(r)(1 − α) ≈ l(r)(1 − α) + N σ (1 − α). It is to be noted that for large N , to achieve a false alarm probability of α, the decision maker requires Nα < N samples (the mean of the log–likelihood ratio, LLR of received samples, after change, is the KL divergence between pdfs f1 and f0 , given by I(f1 , f0 ) > 0. Hence, the posterior probability, which is a function of LLR, increases with the the number of received samples. Thus, to cross a threshold of γ(α), we need Nα samples). Thus, for large N , in the NADM procedure, the detection delay is approximately l(r)(1 − α) + Nσα (1 − α), where Nα /σ is the mean network–delay to transport Nα samples. Thus, for large N , the difference in the mean detection delay between NODM and NADM procedures is approximately 1−α σ (N − Nα ). Note that Nα depends only on α and hence the quantity 1−α (N − N ) increases with N . α σ This behaviour is in agreement with Figure 17. Also, as N × r = 1/3, we expect the network delay to be very large (as 1/3 is close to σ = 0.3636) and hence having a single node is optimal which is also evident from Figure 17. It is also possible to find an example where the optimal number of nodes is greater than 1. For example this occurs in the above setting for N × r = 0.01 (see Figure 18). Note that having N = 10 sensors is optimal for the NADM procedure. The NODM procedure makes the decision only when it receives a batch of N samples corresponding to a sampling instant, whereas NADM procedure makes the decision ACM Transactions on Sensor Networks, Vol. V, No. N, Month 20YY. · Mean Detection Delay (slots) 2800 29 NODM: Simulation NADM: Simulation 2600 2400 2200 2000 1800 1600 1400 1200 1000 800 0 5 10 15 20 25 30 35 40 Number of Sensors, N Fig. 18. Mean detection delay for N × r = 0.01 is plotted against the the number of nodes N . This was obtained with ρ = 0, p = 0.0005, α = 0.01 and with the sensor observations being N (0, 1) and N (1, 1), before and after the change respectively. at every time slot irrespective of whether it receives a sample in that time slot or not. Thus, the Bayesian update that NADM does at every time slot makes it stop earlier than NODM. 7. CONCLUSIONS In this work we have considered the problem of minimizing the mean detection delay in an event detection on a small extent ad hoc wireless sensor network. We provide two ways of processing samples in the fusion center: i) Network Oblivious (NODM) processing, and ii) Network Aware (NADM) processing. We show that in the NODM processing, under periodic sampling, the detection delay decouples into decision and network delays. An important implication of this is that an optimal sequential change detection algorithm can be used in the decision device independently of the random access network. We also formulate and solve the change detection problem in the NADM setting in which case the optimal decision maker needs to use the network state in its optimal stopping rule. Also, we study the network delay involved in this problem and show that it is important to operate at a particular sampling rate to achieve the minimum detection delay. Appendix – I Proof: (Theorem 1)   h i  e − Te + Te − K + K − T )I e e −T I e min E U = min E ( U {T >T } {T >T } Πα Πα r r   h   i K e − Te)I e + E = min E (U − T I {T >T } {Te>T } Πα r   i 1 h e (25) − K I{Te>T } + E K r ACM Transactions on Sensor Networks, Vol. V, No. N, Month 20YY. · 30 Note that in Eqn. 25, the first term is the queueing delay, the second term is the coarse sampling delay and the third term is the decision delay (all delays being in slots). Consider the first term, i i h h e − Te)I e E (U e − tK e )I{Te>T } {T >T } = E (UK   X e = b, Db = x x · I b = P T = j, K { >j} r j≥0,b≥0,x≥0 X =   e = b P (Db = x) x · I b P T = j, K { >j} r j≥0,b≥0,x≥0 where we have used the facts that (i) the decision process is based on only what the packets carry and not on their arrival time etc, and (ii) the condition that sampling is done periodically at a known rate r. Assuming the queueing system to be stationary, the above can be written as   i h  X X  e − Te)I e e =b I b   = E (U P (D = x) x P T = j, K { >j} {T >T } r x≥0   = E[D] P Te > T . j,b Note that E[D] is a function of the sampling rate r, and does not depend on the detection policy. Consider the second term of Eqn. 25,       K K E − T I{Te≥T } = E − T I{K≥K} e r r       K K +E − T I{K≥K,S − T I{K≥K,S = E e e 0 =1} 0 =0} r r For S0 = 1, we have T = 0 and K = 0. Hence,       K K E − T I{Te≥T } = 0 + E0 − T I{K≥K} e r r where E0 [·] denote the expectation and P0 (·) the probability law, when the initial state is S0 = 0. Now, E0  K −T r  I{K≥K} e  = = ∞ X ∞ X b/r X b=1 e b=b t=(b−1)/r+1 ∞ ∞  XX b=1 e b=b  ·     e = eb · b − t P0 T = t, K = b, K r  e = eb P0 K = b, K b/r X t=(b−1)/r+1 P0     b e = eb · T = t | K = b, K − t  (26) r ACM Transactions on Sensor Networks, Vol. V, No. N, Month 20YY. · 31 e is independent of T given K. Hence, We note that K    ∞ X ∞   X K e e E0 K = b, K = b = P − T I{K≥K} 0 e r b=1 e b=b · h 1/r−1 X y=0 We have P0 (T = t | K = b) = (  i b y · P0 T = − y | K = b r (1−ρ)(1−p)t−1 p (1−ρ)(1−pr )b−1 pr , 0, for t s.t. b = ⌈t · r⌉ otherwise. Hence, for 0 ≤ y ≤ 1/r − 1,   b (1 − p)b/r−y−1 p P0 T = − y | K = b = r (1 − pr )b−1 pr But, (1 − pr ) = (1 − p)1/r . Hence,   b (1 − p)b/r−y−1 p P0 T = − y | K = b = r (1 − pr )b−1 pr = (1 − p)1/r−y−1 p 1 − (1 − p)1/r It can be shown that 1/r−1 X y=0 y· (1 − p)1/r−y−1 p 1 − = r 1 − (1 − p)1/r   1 1 (1 − pr ) − p rpr =: l(r) Therefore, Eqn. 26 can be written as          K e ≥ K = l(r) · P K e ≥K −ρ = l(r) · P0 K E0 − T I{K≥K} e r     e <K −ρ = l(r) · 1 − P K Finally, we have i h e − T )I e min E (U {T ≥T } Πα  i     1 h  + e e e + l(r)P0 T ≥ T + E (K − K) = min d(r) 1 − P T < T Πα r     i h 1 + e e − ρ · l(r) + E (K − K) = min (d(r) + l(r)) 1 − P T < T Πα r    is Note that, in the above equation, the first term (d(r) + l(r)) 1 − P Te < T ACM Transactions on Sensor Networks, Vol. V, No. N, Month 20YY. 32 ·   minimum when P Te < T = α. It follows that i h e − T )I e min E (U {T ≥T } Πα > (d(r) + l(r)) (1 − α) − ρ · l(r) + i h 1 e − K)+ min E (K r Πα i   h  e − K)+ achieves 1 − P Te < T = Also, since the optimal policy for the problem minΠα E (K α, we also have (d(r) + l(r)) (1 − α) − ρ · l(r) + h i h i 1 e − K)+ > min E (U e − T )I e min E (K {T ≥T } Πα r Πα It follows that i h h i 1 e − K)+ e − T )I e min E ( K min E (U = (d(r) + l(r)) (1 − α) − ρ · l(r) + {T ≥T } Πα r Πα We need 1 − α > ρ or α < 1 − ρ. If α > 1 − ρ, the optimal stopping isiat t = 0. h e This will yield the desired probability of false alarm and E (U − T )I{Te≥T } = 0. Appendix – II Proof: (Theorem 2) The necessity of N r < σ is clear. The sufficiency proof goes as follows. Consider the FJQ-GPS system with every queue always containing a single dummy packet that is served at low priority. Let us call this the saturated FJQ-GPS system. When a queue becomes empty, the low priority dummy packet contends for service. If it receives service, then it immediately reappears and continues to contend for service. If, while a dummy packet is in service, a regular packet arrives, then the service of the dummy packet is preempted and the regular packet starts contending. It follows that the service rate applied to every queue (i.e., those with regular packets or those with dummy packets) is always σ/N . Now, consider a virtual service process of rate σ. In each slot, a service occurs with probability σ and the service is applied to any one of the queues with equal probability. Equivalently each queue is served by an independent Bernoulli process of rate σ/N . Considering only the services to the regular packets at each queue, we have a GI/M/1 queue (here GI refers to a General distribution with Independent arrivals, M refers to a Markovian service process and 1 refers to one server). Hence, the system has proper stationary delay, iff r < σ/N . Also, it can be seen that the delays in the above described system (with dummy packets when a queue is empty) upper bound those in the original FJQ-GPS system. Hence, the result follows. Appendix – III Distribution of state noise N  

Let q = [λ, b, δ, w, r]. Note that P Mk = mQk = q, Θk = θ = P Mk = mQk = q ACM Transactions on Sensor Networks, Vol. V, No. N, Month 20YY. · 33 and is given by 
P Mk = 0Qk = q = 
P Mk = mQk = q =   1 1−σ if φN (q) = 0 if φN (q) > 0 0 if φN (q) = 0 if φL(m) (q) > 0, m = 1, 2, 3, · · · , N. σ φN (q) where φN (q) and φL(m) (q) are obtained from Eqns. 9 and 8.  

The distribution function, P Ok = oQk = q, Θk = θ = P Ok = oQk = q, Θk = θ is given by  
1−p P Ok = oQk = q, Θk = 0 = p  0  
1 P Ok = oQk = q, Θk = 1 = 0 if o = 0 if o = 1, otherwise. if o = 0 otherwise. Appendix – IV Proof of Lemma–1 Let q = [λ, b, δ, w, r]. From Eqn. 17, 
Πk := P T 6 k Ik−1 , Qk = q, Yk  

= P T 6 k − δIk−1 , Qk = q, Yk + P k − δ < T 6 k Ik−1 , Qk = q, Yk   


= P T 6 k − δIk−1 , Qk = q, Yk + P T > k − δIk−1 , Qk = q, Yk · P T 6 k T > k − δ, Ik−1 , Qk = q, Yk , 
= Ψk + (1 − Ψk ) · P T 6 k T > k − δ, Ik−1 , Qk = q, Yk , 
P (k − δ < T 6 k) P Ik−1 , Qk = q, Yk k − δ < T 6 k 
= Ψk + (1 − Ψk ) · P (T > k − δ) P Ik−1 , Qk = q, Yk T > k − δ P (k − δ < T 6 k) P (T > k − δ)   = Ψk + (1 − Ψk ) 1 − (1 − p)δ = Ψk + (1 − Ψk ) · ACM Transactions on Sensor Networks, Vol. V, No. N, Month 20YY. (27) (28) 34 · Eqn. 27 is justified as follows. Note that 
P Ik−1 , Qk = q, Yk k − δ < T 6 k    (i) (i) = P Q[0:k−1] , Qk = q, X[1:Bk −1] , {XBk : Rk = 1}, u[0:k−1]k − δ < T 6 k 
= P Q[0:k−1] , Qk = qk − δ < T 6 k    (i) (i) ·P X[1:Bk −1] , {XBk : Rk = 1}k − δ < T 6 k, Q[0:k−1] , Qk = q    (i) (i) ·P u[0:k−1] k − δ < T 6 k, Q[0:k−1] , Qk = q, X[1:Bk −1] , {XBk : Rk = 1}   
(i) (i) = P Q[0:k−1] , Qk = q · P X[1:Bk −1] , {XBk : Rk = 1}k − δ < T, Q[0:k−1] , Qk = q    (i) (i) ·P u[0:k−1] Q[0:k−1] , Qk = q, X[1:Bk −1] , {XBk : Rk = 1}    
(i) (i) = P Q[0:k−1] , Qk = qT > k − δ · P X[1:Bk −1] , {XBk : Rk = 1}T > k − δ, Q[0:k−1] , Qk = q    (i) (i) ·P u[0:k−1] T > k − δ, Q[0:k−1] , Qk = q, X[1:Bk −1] , {XBk : Rk = 1} 
= P Ik−1 , Qk = q, Yk T > k − δ . We use the following facts in the above justification: i) the evolution of the queueing system Qk is independent of the change point T , ii) whenever T > k − δ, the (i) distribution of any sample Xh , h 6 Bk is f0 , and iii) the control uk = µ̃(Ik ). Appendix – V Recursive computation of Πk At time k, based on the index of the node that successfully transmits a packet Mk , the set of all sample paths Ω can be partitioned based on the following events, n o (j) E1,k := ω : Mk (ω) = 0 or Mk (ω) = j > 0, Rk (ω) = 1 ) ( N X (i) (j) Rk (ω) < N − 1 E2,k := ω : Mk (ω) = j > 0, Rk (ω) = 0, E3,k := ( i=1 ω : Mk (ω) = j > (j) 0, Rk (ω) = 0, N X (i) Rk (ω) i=1 ) =N −1 , i.e., Ω = E1,k ∪ E2,k ∪ E3,k . We note that the above events can also be described by using Qk and Qk+1 in the following manner E1,k = {ω : Wk+1 (ω) = Wk (ω), Rk+1 (ω) = Rk (ω)} [ {ω : Wk+1 (ω) = Wk (ω) + ej , Rk+1 (ω) = Rk (ω)} E2,k = {ω : Wk+1 (ω) = Wk (ω), Rk+1 (ω) = Rk (ω) + ej } ) ( N X (i) (i) (i) (i) + Rk (ω) = N − 1, ∀i, Wk+1 (ω) = (Wk (ω) − 1) , Rk+1 (ω) = 1{W (i) >0} . E3,k = ω : i=1 k Here, the events E1,k and E2,k represent the case Bk+1 = Bk , and the event E3,k represents the case Bk+1 = Bk + 1 (i.e., only if the event E3,k occurs then the ACM Transactions on Sensor Networks, Vol. V, No. N, Month 20YY. · 35 batch index is incremented). We are interested in obtaining Πk+1 from [Qk , Πk ] and Zk+1 . We show that at time k + 1, the statistic Ψk+1 (after having observed Zk+1 ) can be computed in a recursive manner using Ψk and Qk . Using Lemma 1 (using Eqn. 19) we compute Πk+1 from Ψk+1 .   e k+1 = 1 | Ik+1 Ψk+1 = P Θ = = 3   X e k+1 = 1, Ec,k | Ik+1 P Θ c=1 3 X c=1   e k+1 = 1 | Ec,k , Ik+1 1E P Θ c,k (∵ Ec,k is Ik+1 measurable) (j) • Case Mk = 0 or Mk = j > 0, Rk = 1: Πk+1 = P (Θk+1 = 1 | E1,k , Ik+1 ) = P (Θk+1 = 1 | E1,k , Ik , Qk+1 = q′ ) P (Θk+1 = 1 | E1,k , Ik ) · fQk+1 |Θk+1 ,E1,k ,Ik (q′ |1, E1,k , Ik ) =  f (q′ |E1,k , Ik )  Qk+1 = P (Θk+1 = 1 | E1,k , Ik ) (by Bayes rule) E1,k ,Ik (Qk+1 is independent of Θk+1 ) = P (Θk = 0, Θk+1 = 1 | Ik ) + P (Θk = 1, Θk+1 = 1 | Ik ) = (1 − Πk )p + Πk PN (i) (j) (j) • Case Mk = j > 0, Rk = 0, i=1 Rk < N − 1: In this case, the sample XBk is successfully transmitted and is passed on to the decision maker. The decision maker receives just this sample, and computes Πk+1 . We compute Ψk+1 from Ψk and then we use Lemma 1 (using Eqn. 19) to compute Πk+1 from Ψk+1 . Ψk+1   e k+1 = 1 | E2,k , Ik+1 = P Θ   e k+1 = 1 | E2,k , Ik , [Qk+1 , Yk+1 ] = [q′ , y] = P Θ   e k = 0, Θ e k+1 = 1 | E2,k , Ik , [Qk+1 , Yk+1 ] = [q′ , y] = P Θ   e k = 1, Θ e k+1 = 1 | E2,k , Ik , [Qk+1 , Yk+1 ] = [q′ , y] +P Θ Since, we consider the case when the fusion center received a sample at time e k+1 = Θk+1−∆ k + 1 and Bk+1 = Bk , ∆k+1 = ∆k + 1 and hence, the state Θ = k+1 e Θk−∆k = Θk . Thus, in this case, Ψk+1 can be written as ACM Transactions on Sensor Networks, Vol. V, No. N, Month 20YY. · 36 Ψk+1   e k = 1, Θ e k+1 = 1 | E2,k , Ik , [Qk+1 , Yk+1 ] = [q′ , y] = P Θ     e k = 1, Θ e k+1 = 1 | E2,k , Ik · P Qk+1 = q′ | Θ e k = 1, Θ e k+1 = 1, E2,k , Ik P Θ (a) = P(Qk+1 = q′ |E2,k , Ik ) · fYk+1 |E2,k ,Ik ,Qk+1 (y|E2,k , Ik , q′ ) (b) = (y | 1, 1, E2,k , q′ , Ik ) ·fY e e k+1 |Θk ,Θk+1 ,E2,k ,Ik ,Qk+1   e k = 1, Θ e k+1 = 1 | E2,k , Ik · P(Qk+1 = q′ |E2,k , Ik ) · f P Θ e Yk+1 |Θ k q′ |E2,k , Ik ) (y | 1) , q′ ) · fYk+1 |E2,k ,Ik ,Qk+1 (y|E2,k , Ik   e k = 1, Θ e k+1 = 1 | E2,k , Ik · f1 (y) P Θ (c)     = e k = 0 | E2,k , Ik , Qk+1 · f e P Θ e (y|0) + P Θk = 1 | E2,k , Ik , Qk+1 · fY Y |Θ P(Qk+1 = k+1 k e k+1 |Θk (y|1) Ψk f1 (y) = (1 − Ψk )f0 (y) + Ψk f1 (y) (d) We explain the steps (a), (b), (c), (d) below. (a) By Bayes rule, for events A, B, C, D, E, F , we have P (AB | CDEF ) = P (AB | CD) P (E | ABCD) P (F | ABCDE) P (E | CD) P (F | CDE) e k, Θ e k+1 . Also, given Θ e k , Yk+1 is independent of (b) Qk+1 is independent of Θ e Θk+1 , E2,k , Ik , Qk+1 (c) For any events A, B, and a continuous random variable Y , the conditional density function fY |A (y|A) = P (B | A) fY |AB (y|AB)+P (B c | A) fY |AB c (y|AB c ). e k , Yk+1 is independent of E2,k , Ik , Qk+1 Also, given Θ e k is independent (d) E2,k is [Ik , Qk+1 ] measurable, and hence, given [Ik , Qk+1 ], Θ of E2,k . PN (i) (j) • Case Mk = j > 0, Rk = 0, i=1 Rk = N − 1: In this case, at time k + 1, the (j) decision maker receives the last sample of batch Bk , XBk (that is successfully transmitted during slot k) and the samples of batch Bk + 1, if any, that are queued in the sequencer buffer. We compute Ψk+1 from Ψk and then we use Lemma 1 (using Eqn. 19) to compute Πk+1 from Ψk+1 . In this case, the decision PN maker receives n := i=1 1{W (i) >0} samples of batch Bk + 1. Also, note that n k is Ik measurable.   e k+1 = 1 | E3,k , Ik+1 Ψk+1 = P Θ   e k+1 = 1 | E3,k , Ik , [Qk+1 , Yk+1 ] = [q′ , y] = P Θ   e k = 0, Θ e k+1 = 1 | E3,k , Ik , [Qk+1 , Yk+1 ] = [q′ , y] = P Θ   e k = 1, Θ e k+1 = 1 | E3,k , Ik , [Qk+1 , Yk+1 ] = [q′ , y] +P Θ Since, we consider the case Bk+1 = Bk + 1, ∆k+1 = ∆k + 1 − 1/r and hence, the e k+1 = Θk+1−∆ state Θ = Θk−∆k +1/r . k+1 ACM Transactions on Sensor Networks, Vol. V, No. N, Month 20YY. · 37 Let y = [y0 , y1 , · · · , yn ]. Consider (a) = (b) =   eΘ e k = θ, e k+1 = 1 | E3,k , Ik , [Qk+1 , Yk+1 ] = [q′ , y] P Θ     eΘ eΘ e k = θ, e k+1 = 1 | E3,k , Ik · P Qk+1 = q′ | Θ e k = θ, e k+1 = 1, E3,k , Ik P Θ P(Qk+1 = q′ |E3,k , Ik ) · fYk+1 |E3,k ,Ik ,Qk+1 (y|E3,k , Ik , q′ ) e 1, E3,k , q′ , Ik ) ·fY (y | θ, e e k+1 |Θk ,Θk+1 ,E3,k ,Ik ,Qk+1   eΘ e k = θ, e k+1 = 1 | E3,k , Ik · P(Qk+1 = q′ |E3,k , Ik ) · f e(y0 ) Qn f1 (yi ) P Θ i=1 θ P(Qk+1 = q′ |E3,k , Ik ) · fYk+1 |E3,k ,Ik ,Qk+1 (y|E3,k , Ik , q′ )     e E3,k , Ik · f e(y0 ) Qn f1 (yi ) e k+1 = 1 | Θ e k = θ, e k = θe | E3,k , Ik · P Θ P Θ i=1 θ (c)   . = P P1 1 ′ ′′ e e e e (y|θe′ , θe′′ , E3,k , Ik , q′ ) e ,Θ e e′ =0 e′′ =0 P Θk = θ , Θk+1 = θ , | E3,k , Ik , Qk+1 · fY |Θ E ,I ,Q θ θ k+1 k k+1 3,k k k+1 We explain the steps (a), (b), (c) below. (a) By Bayes rule, for events A, B, C, D, E, F , we have P (AB | CDEF ) = P (AB | CD) P (E | ABCD) P (F | ABCDE) P (E | CD) P (F | CDE) e k, Θ e k+1 . Also, given Θ e k , Yk+1,0 is independent of (b) Qk+1 is independent of Θ e e e k , E3,k , Ik , Qk+1 . Θk+1 , E3,k , Ik , Qk+1 , and given Θk+1 , Yk+1,i is independent of Θ It is to be noted that given the state of nature, the sensor measurements Yk+1,0 , Yk+1,1 , · · · , Yk+1,n are conditionally independent. (c) For any events A, B, and a continuous random variable Y , the conditional density function fY |A (y|A) = P (B | A) fY |AB (y|AB)+P (B c | A) fY |AB c (y|AB c ). e k , Yk+1 is independent of E3,k , Ik , Qk+1 Also, given Θ It is to be noted that the event E3,k is [Ik , Qk+1 ] measurable, and hence, given e k is independent of E3,k . Thus, in this case, [Ik , Qk+1 ], Θ Ψk+1 = Q Qn (1 − Ψk )pr f0 (y0 ) n i=1 f1 (yi ) + Ψk f1 (y0 ) i=1 f1 (yi ) Qn Qn Q . (1 − Ψk )(1 − pr )f0 (y0 ) i=1 f0 (yi ) + (1 − Ψk )pr f0 (y0 ) i=1 f1 (yi ) + Ψk f1 (y0 ) n i=1 f1 (yi ) Thus, using Lemma 1 (using Eqn. 19), we have Πk+1 = Ψk+1 + (1 − Ψk+1 )(1 − (1 − p)∆k+1 ) =: φΨ (Ψk , Zk+1 ) + (1 − φΨ (Ψk , Zk+1 )) (1 − (1 − p)∆k+1 )   Πk − (1 − (1 − p)∆k ) , Z = φΨ k+1 (1 − p)∆k    Πk − (1 − (1 − p)∆k ) (1 − (1 − p)∆k+1 ) , Zk+1 + 1 − φΨ (1 − p)∆k =: φΠ ([Qk , Πk ], Zk+1 ) . Appendix – VI Structure of τ ∗ We use the following Lemma to show that J ∗ (q, π) is concave in π. ACM Transactions on Sensor Networks, Vol. V, No. N, Month 20YY. 38 · Lemma 2 If f : [0, 1] → R is concave, then the function h : [0, 1] → R defined by    y · φ2 (x) + (1 − y)pr · φ1 (x) h(y) = Eφ(x) f y · φ2 (x) + (1 − y)pr · φ1 (x) + (1 − y)(1 − pr ) · φ0 (x) is concave for each x, where φ(x) = y·φ2 (x)+(1−y)pr ·φ1 (x)+(1−y)(1−pr )·φ0 (x), 0 < pr < 1, and φ0 (x), φ1 (x), and φ2 (x) are pdfs on X. Proof. See Appendix – I of [Premkumar and Kumar 2008]. Note that in the finite H–horizon (truncated version of Eqn. 21), we note from value H iteration that the cost–to–go function, for a given q, JH ([q, π]) = 1 − π is concave in π. Hence, by Lemma 2, we see that for any given q, the cost–to–go functions H H JH−1 ([q, π]), JH−2 ([q, π]), · · · , J0H ([q, π]) are concave in π. Hence for 0 ≤ λ ≤ 1, J ∗ ([q, π]) = lim J0H ([q, π]) H→∞   ∗ J ([q, λπ1 + (1 − λ)π2 ]) = lim J0H [q, λπ1 + (1 − λ)π2 ] H→∞ ≥ lim λJ0H ([q, π1 ]) + lim (1 − λ)J0H ([q, π2 ]) H→∞ ∗ H→∞ ∗ = λJ ([q, π1 ]) + (1 − λ)J ([q, π2 ]) It follows that for any given q, J ∗ ([q, π]) is concave in π. Define the map ξ : Q×[0, 1] → R+ as ξ([q, π]) Q×[0, 1] →   := 1−π and the map κ :   ∗ R+ , as κ([q, π]) := c·π+AJ ∗ ([q, π]) = c·π+E J ([Qk+1 , φΠ (νk , Zk+1 )])  νk = [q, π] . Note that ξ([q, 1]) = 0, κ([q, 1]) = c, ξ([q, 0]) = 1 and     κ([q, 0]) = E J ∗ ([Qk+1 , φΠ (νk , Zk+1 )])  ν = [q, 0] k      (2) ∗ = E J ([φQ (Qk , Mk ), φΠ (νk , Zk+1 )])  ν = [q, 0] k  = N X m=0 (4) 6 N X m=0 = N X m=0 (6) 6 N X m=0         ∗  E J ([φQ (q, m), φΠ (νk , Zk+1 )]) Mk = m, νk = [q, 0] P Mk = m ν = [q, 0] k  J ∗           φQ (q, m), E φΠ (νk , Zk+1 )Mk = m, νk = [q, 0] P Mk = mνk = [q, 0]      J ∗ ([φQ (q, m), p) P Mk = m ν = [q, 0]  k     (1 − p) · P Mk = m ν = [q, 0]  k = 1−p < 1 where in the above derivation, we use the evolution of Qk in step 2, the Jensen’s inequality (as for any given q, J ∗ (q, π) is concave in π) in step 4, and the inequality J ∗ (q, π) 6 1 − π in step 6. ACM Transactions on Sensor Networks, Vol. V, No. N, Month 20YY. · 39 Note that κ([q, 1]) − ξ([q, 1]) > 0 and κ([q, 0]) − ξ([q, 0]) < 0. Also, for a fixed q, the function κ([q, π]) − ξ([q, π]) is concave in π. Hence, by the intermediate value theorem, for a fixed q, there exists γ(q) ∈ [0, 1] such that κ([q, γ]) = ξ([q, γ]). This γ is unique as κ([q, π]) = ξ([q, π]) for at most two values of π. If in the interval [0, 1], there are two distinct values of π for which κ([q, π]) = ξ([q, π]), then the signs of κ([q, 0]) − ξ([q, 0]) and κ([q, 1]) − ξ([q, 1]) should be the same. Hence, τ ∗ = inf {k : Πk > γ(Qk )} where the threshold γ(q) is given by c · γ(q) + AJ ∗ ([q, γ(q)]) = 1 − γ(q). REFERENCES Aldosari, S. A. and Moura, J. M. F. 2004. Detection in decentralized sensor networks. In Proceedings of ICASSP. II:277–280. Baccelli, F. and Makowski, A. 1990. Synchronization in queueing systems. In Stochastic Analysis of Computer and Communication Systems, H. Takagi, Ed. North–Holland, 57–131. Bertsekas, D. P. 2000a. Dynamic Programming and Optimal Control , Second ed. Vol. I. Athena Scientific. Bertsekas, D. P. 2000b. Dynamic Programming and Optimal Control , Second ed. Vol. II. Athena Scientific. Honeywell Inc . http://hpsweb.honeywell.com/Cultures/en-US/Products/Wireless/ SecondGenerationWireless/default.htm. ISA . http://www.isa.org/Content/NavigationMenu/Technical_Information/ASCI/ISA100_ Wireless_Compliance_Institute/ISA100_Wireless_Compliance_Institute.htm. Kumar, A., Manjunath, D., and Kuri, J. 2004. Communication Networking: An Analytical Approach. Morgan-Kaufmann (an imprint of Elsevier), San Francisco. Kumar, A., Manjunath, D., and Kuri, J. 2008. Wireless Networking. Morgan-Kaufmann (an imprint of Elsevier), San Francisco. Niu, R. and Varshney, P. K. 2005. Distributed detection and fusion in a large wireless sensor network of random size. EURASIP Journal on Wireless Communications and Networking 4, 7, 462–472. Premkumar, K. and Kumar, A. 2008. Optimal sleep-wake scheduling for quickest intrusion detection using wireless sensor networks. In Proc. IEEE Infocom. AZ, USA. Premkumar, K., Kumar, A., and Kuri, J. 2009. Distributed detection and localization of events in large ad hoc wireless sensor networks. In Proc. 47th Annual Allerton Conference on Communication, Control, and Computing. IL, USA. Rajasegarar, S., Leckie, C., and Palaniswami, M. 2008. Anomaly detection in wireless sensor networks. IEEE Wireless Communications 15, 4 (August), 34–40. Shiryaev, A. N. 1978. Optimal Stopping Rules. Springer, New York. Singh, C. K., Kumar, A., and Ameer, P. M. 2008. Performance evaluation of an IEEE 802.15.4 sensor network with star topology. Wireless Networks 14, 4 (August), 543–568. Tartakovsky, A. G. and Veeravalli, V. V. 2005. General asymptotic bayesian theory of quickest change detection. SIAM Theory of Probability and its Applications 49, 3, 458–497. Tenny, R. R. and Sandell, N. R. 1981. Detection with distributed sensors. IEEE Transactions on Aerospace and Electronic Systems 17, 501–510. Veeravalli, V. V. 2001. Decentralized quickest change detection. IEEE Transactions on Information theory 47, 4 (May), 1657–1665. Received March 2009; revised December 2009 and June 2010; accepted Month Year ACM Transactions on Sensor Networks, Vol. V, No. N, Month 20YY.