Dr Weihua Yang

Dr Yandong Bai

Dr Weihua He

Dr Qiang Sun

Dr Kai Yang

Dr Guangyu Li

Ms Mr Chuan Duzhong Xu

Zhang

First of all, I would like to thank my advisor, Dr. Lin Chen, who illuminates my way to science. Thanks to his open-door policy, I can always discuss

Introduction 1.1 Background and Motivation

Recent years have witnessed an unprecedented development of the radio frequency identification (RFID) technology due to its low cost and non-line-of-sight communication pattern which overcomes the drawbacks of traditional barcode technology [1]. An RFID system typically consists of one or several RFID readers and a large number of RFID tags. Specifically, an RFID reader is a device equipped with a dedicated power source and an antenna and can collect and process the information of tags within its coverage area. An RFID tag, on the other hand, is a low-cost microchip labeled with a unique serial number (ID) to identify an object and can receive and transmit the radio signals via the wireless channel. More specially, the tags are generally classified into two categories: passive tags and active tags. The passive tags are energized by the radio wave of the reader, whereas, the active tags have power sources and relatively long communication range. Moreover, the communication between readers and tags follows frame slotted Aloha (FSA) mechanism which has been standardized in one of the most popular industrial standards, EPCGlobal Class-1 Generation-2 (C1G2) RFID standard [2].

As a promising low-cost technology, RFID is widely utilized in various applications ranging from inventory control to supply chain management and logistics. In most, if not all, RFID applications, tag counting and monitoring are perhaps one of the most fundamental component. Although simple to state and intuitively understandable, designing efficient tag counting and monitoring algorithms require non-trivial efforts to solve, especially in large-scale RFID systems, due to the following particular challenges in RFID systems.

• Large number of tags. An RFID system may consist of a large number of tags, such as a warehouse storing thousands of goods for retailers. Any effective algorithm designed for these RFID systems needs to scale elegantly.

• Limited computing resource at tags. The quest of compatible size and low energy consumption significantly limits the computing and processing capability of individual tags in RFID systems, especially for lightweight passive tags.

• Unreliable wireless links. Wireless links are notoriously unreliable and error-prone. Hence, algorithms To lay the theoretical foundations for the design and optimization of tag counting and monitoring algorithms, we start by investigating the stability of FSA, which is of fundamental importance both on the theoretical characterisation of FSA performance and its effective operation in practical systems. To study the stability of FSA, the effort should be devoted to answering the following natural and crucial questions:

• Under what condition(s) is FSA stable?

• When is the stability region maximised?

• How does FSA behave in the instability region? Surprisingly, the above fundamental questions have not been explored in the literature. To fill the void in the study of FSA stability, in Chapter 2, we investigate these questions by technically modeling the FSA system backlog as a Markov chain with its states being backlog size at the beginning of each frame. The main objective is to analyze the ergodicity of the Markov chain and demonstrate its properties in different regions, particularly the instability region. Specifically, by employing drift analysis, we obtain the closed-form conditions for the stability of FSA and show that the stability region is maximised when the frame length equals the number of packets to be sent in the single packet reception model and the upper bound of stability region is maximised when the ratio of the number of packets to be sent to frame length equals in an order of magnitude the maximum multipacket reception capacity in the multipacket reception model, which answers the first two questions. Furthermore, to characterise system behavior in the instability region, we mathematically demonstrate the existence of transience of the backlog Markov chain, which provides the answer to the third question.

Chapter 3: From Static to Dynamic Tag Population Estimation: An Extended Kalman Filter Perspective

In this chapter, we focus on the problem of tag counting, or tag population estimation, which has recently attracted significant research attention due to its paramount importance on a variety of RFID applications. However, most, if not all, of existing estimation mechanisms are proposed for the static case where tag population remains constant during the estimation process, thus leaving the more challenging dynamic case unaddressed, despite the fundamental importance of the latter case on both theoretical analysis and practical application.

Motivated by the above argument, we design a generic framework of stable and accurate tag population estimation schemes based on Kalman filter for both static and dynamic RFID systems. Our main contributions are two-fold. Firstly, we model the dynamics of RFID systems as discrete stochastic processes and leverage the techniques in extended Kalman filter (EKF) and cumulative sum control chart (CUSUM) to estimate tag population for both static and dynamic systems. Secondly, By employing Lyapunov drift analysis, we mathematically characterise the performance of the proposed framework in terms of estimation accuracy and convergence speed by deriving the closed-form conditions on the design parameters under which our scheme can stabilise around the real population size with bounded relative estimation error that tends to zero with exponential convergence rate.

Chapter 4: Finding Needles in a Haystack: Missing Tag Detection in the Presence of Unexpected Tags RFID technology has been widely used in missing tag detection to reduce and avoid inventory shrinkage. In this application, promptly finding out the missing event is of paramount importance. In Chapter 4 and 5, we further address the problem of missing tag detection. Different to existing works in this field, we focus on two unexplored while fundamentally important scenarios, missing tag detection in the presence of unexpected tags (Chapter 4) and in multiple-group multiple-region RFID systems (Chapter 5).

In the first scenario, existing missing tag detection protocols cannot efficiently handle the presence of a large number of unexpected tags whose IDs are not known to the reader, which shackles the time efficiency. To deal with the problem of detecting missing tags in the presence of unexpected tags, we devise a two-phase Bloom filter-based missing tag detection protocol (BMTD). The proposed BMTD exploits Bloom filter in sequence to first deactivate the unexpected tags and then test the membership of the expected tags, thus dampening the interference from the unexpected tags and considerably reducing the detection time. Moreover, the theoretical analysis of the protocol parameters is performed to minimize the detection time of the proposed BMTD and achieve the required reliability simultaneously. Extensive experiments are then conducted to evaluate the performance of the proposed BMTD. The results demonstrate that the proposed BMTD significantly outperforms the state-of-the-art solutions.

Chapter 5: On Missing Tag Detection in Multiple-group Multiple-region RFID Systems

In Chapter 5, we formulate and study a new missing tag detection problem, arising in multiple-group multipleregion RFID systems, where a mobile reader needs to detect whether there is any missing event for each group of tags. The problem we tackle is to devise missing tag detection protocols with minimum execution time while guaranteeing the detection reliability requirement for each group. By leveraging the technique of Bloom filter, we develop a suite of three missing tag detection protocols, each decreasing the execution time compared to its predecessor by incorporating an improved version of the Bloom filter design and parameter tuning. By sequentially analysing the developed protocols, we gradually iron out an optimum detection protocol that works in practice.

Finally, Chapter 6 concludes the thesis with the summary of the overall results and the perspective for the future research.

Part of our research work presented in this thesis is published or in submission in various journals and conferences. Specifically, our work on the stability analysis of FSA was presented in the 23rd International

Context and Motivation

Since the introduction of Aloha protocol in 1970 [3], a variety of such protocols have been proposed to improve its performance, such as Slotted Aloha (SA) [4] and Frame Slotted Aloha (FSA) [START_REF] Okada | Analysis and application of framed ALOHA channel in satellite packet switching networks-fadra method[END_REF]. SA is a well known random access scheme where the time of the channel is divided into identical slots of duration equal to the packet transmission time and the users contend to access the server with a predefined slot-access probability. As a variant of SA, FSA divides time-slots into frames and a user is allowed to transmit only a single packet per frame in a randomly chosen time-slot.

Due to their effectiveness to tackle collisions in wireless networks, SA-and FSA-based protocols have been applied extensively to various networked systems ranging from the satellite networks [START_REF] Noh | Anc-aloha: Analog network coding ALOHA for satellite networks[END_REF], wireless LANs [START_REF] Vasudevan | Neighbor discovery in wireless networks and the coupon collector's problem[END_REF][START_REF] Zeng | Neighbor discovery in wireless networks with multipacket reception[END_REF] to the emerging Machine-to-Machine (M2M) networks [START_REF] Wu | Fasa: Accelerated S-ALOHA using access history for event-driven M2M communications[END_REF][START_REF] Vázquez | Energy and delay analysis of contention resolution mechanisms for machine-to-machine networks based on low-power wifi[END_REF]. Specifically, in radio frequency identification (RFID) systems, which is our specific interest in this thesis, FSA plays a fundamental role in the identification of tags [START_REF] Zhu | Energy-efficient identification in large-scale rfid systems with handheld reader[END_REF][START_REF] Liu | Efficient unknown tag identification protocols in large-scale rfid systems[END_REF] and is standardized in the EPCGlobal Class-1 Generation-2 (C1G2) RFID standard [2]. In FSAbased protocols, all users with packets transmit in the selected slot of the frame respectively, but only packets experiencing no collisions are successful while the other packets referred to as backlogged packets (or simply backlogs), are retransmitted in the subsequent frames.

Given the paramount importance of the stability for systems operating on top of Aloha-based protocols, a large body of studies have been devoted to stability analysis in a slotted collision channel [START_REF] Bordenave | Asymptotic stability region of slotted aloha[END_REF][START_REF] Kompalli | On the stability of finite queue slotted aloha protocol[END_REF][START_REF] Farhadi | Stability region of a slotted aloha network with k-exponential backoff[END_REF] where a transmission is successful if and only if just a single user transmits in the selected slot, referred to as single packet reception (SPR). Differently with SPR, the emerging multipacket reception (MPR) technologies in wireless networks, such as Code Division Multiple Access (CDMA) and Multiple-Input and Multiple-Output (MIMO), make it possible to receive multiple packets in a time-slot simultaneously, which remarkably boosts system performance at the cost of the system complexity.

More recently, the application of FSA in RFID systems and M2M networks has received considerable research attention. However, very limited work has been done on the stability of FSA despite its fundamental importance both on the theoretical characterisation of FSA performance and its effective operation in practical systems. Motivated by the above observation, we argue that a systematic study on the stability properties of FSA incorporating the MPR capability is called for in order to lay the theoretical foundations for the design and optimization of FSA-based communication systems.

Summary of Contributions

In this chapter, we investigate the stability properties of p-persistent FSA with SPR and MPR capabilities. The main contributions of this chapter are articulated as follows:

• We model the packet transmission process in a frame as the bins and balls problem [START_REF] Johnson | Urn models and their application: an approach to modern discrete probability theory[END_REF] and derive the number of successfully received packets under both SPR and MPR models.

• We formulate a homogeneous Markov chain to characterize the number of the backlogged packets and derive the one-step transition probability with the persistence probability p.

• By employing drift analysis, we obtain the closed-form conditions for the stability of p-persistent FSA and derive conditions maximising the stability regions for both SPR and MPR models.

• To characterise system behavior in the instability region with the persistence probability p, we mathematically demonstrate the existence of transience of the backlog Markov chain.

• We investigate how to achieve the stability condition and give the control algorithm for updating the frame size.

Our work demonstrates that the stability region is maximised when the frame length equals the number of sent packets in the SPR model and the upper bound of stability region is maximised when the ratio of the number of sent packets to frame length equals in an order of magnitude the maximum multipacket reception capacity in the MPR model. In addition, it is also shown that FSA-MPR outperforms FSA-SPR remarkably in terms of the stability region size.

Related Work

Aloha-based protocols are basic schemes for random medium access and are applied extensively in many communication systems. As a central property, the stability of Aloha protocols has received a lot of research attention, which we briefly review in this section.

Stability of slotted Aloha. Tsybakov and Mikhailov [START_REF] Tsybakov | Ergodicity of a slotted ALOHA system[END_REF] initiated the stability analysis of finite-user slotted Aloha. They found sufficient conditions for stability of the queues in the system using the principle of stochastic dominance and derived the stability region for two users explicitly. For the case of more than two users, the inner bounds to the stability region were shown in [START_REF] Rao | On the stability of interacting queues in a multiple-access system[END_REF]. Subsequently, Szpankowski [START_REF] Szpankowski | Stability conditions for some distributed systems: Buffered random access systems[END_REF] found necessary and sufficient conditions for the stability under a fixed transmission probability vector for three-user case. However, the derived conditions are not closed-form, meaning the difficulty on verifying them. In [START_REF] Bordenave | Asymptotic stability region of slotted aloha[END_REF] an approximate stability region was derived for an arbitrary number of users based on the mean-field asymptotics. It was claimed that this approximate stability region is exact under large user population and it is accurate for smallsized networks. The sufficient condition for the stability was further derived to be linear in arrival rates without the requirement on the knowledge of the stationary joint statistics of queue lengths in [START_REF] Kompalli | On the stability of finite queue slotted aloha protocol[END_REF]. Recently, the stability region of SA with K-exponential backoff was derived in [START_REF] Farhadi | Stability region of a slotted aloha network with k-exponential backoff[END_REF] by modeling the network as inter-related quasibirth-death processes. We would like to point out that all the above stability analysis results were derived for the SPR model.

Stability of slotted Aloha with MPR.

The first attempt at analyzing stability properties of SA with MPR was made by Ghez et al. in [START_REF] Ghez | Stability properties of slotted Aloha with multipacket reception capability[END_REF][START_REF] Ghez | Optimal decentralized control in the random access multipacket channel[END_REF] in an infinite-user single-buffer model. They drew a conclusion that the system could be stabilized under the symmetrical MPR model with a non-zero probability that all packets were transmitted successfully. Afterwards, Sant and Sharma [START_REF] Sant | Performance analysis of a slotted-ALOHA protocol on a capture channel with fading[END_REF] studied a special case of the symmetrical MPR model for finite-user with an infinite buffer. They derived sufficient conditions on arrival rate for stability of the system under the stationary ergodic arrival process. Subsequently, the effect of MPR on stability and delay was investigated in [START_REF] Naware | Stability and delay of finite-user slotted ALOHA with multipacket reception[END_REF] and it was shown that stability region undergoes a phase transition and then reaches the maximization. Besides, in [START_REF] Inaltekin | Selfish random access over wireless channels with multipacket reception[END_REF] necessary and sufficient conditions are obtained for a Nash equilibrium strategy for wireless networks with MPR based on noncooperation game theory. More recently, Jeon and

Ephremides [START_REF] Jeon | On the stability of random multiple access with stochastic energy harvesting[END_REF] characterised the exact stability region of SA with stochastic energy harvesting and MPR for a pair of bursty users. Although the work aforementioned analyzed the stability of system without MPR or/and with MPR, they are mostly, if not all, focused on SA protocol, while our focus is FSA with both SPR and MPR.

Performance analysis of FSA.

There exist several studies on the performance of FSA. Wieselthier and Anthony [START_REF] Wieselthier | An exact analysis and performance evaluation of framed ALOHA with capture[END_REF] introduced an combinatorial technique to analyse performance of FSA-MPR for the case of finite users. Schoute [START_REF] Schoute | Dynamic frame length ALOHA[END_REF] investigated dynamic FSA and obtained the expected number of time-slots needed until the backlog becomes zero. Recently, the optimal frame setting for dynamic FSA was proved mathematically in [START_REF] Prodanoff | Optimal frame size analysis for framed slotted Aloha based rfid networks[END_REF] and [START_REF] Barletta | A formal proof of the optimal frame setting for dynamic-frame aloha with known population size[END_REF]. However, these works did not address the stability of FSA, which is of fundamental importance.

In summary, only very limited work has been done on the stability of FSA despite its fundamental importance both on the theoretical characterisation of FSA performance and its effective operation in practical systems. In order to bridge this gap, we devote this chapter to investigating the stability properties of FSA under both SPR and MPR models.

System Model

We introduce our system model which will be used throughout the rest of this chapter.

Physical layer and random access model in FSA

We consider a system of infinite identical users operating on one frequency channel. In one slot, a node can complete a packet transmission. We investigate two physical layer models of practical importance, the models with single packet reception (SPR) and multipacket reception (MPR) capabilities:

• Under the SPR model, a packet suffers a collision if more than one packet is transmitted in the same time-slot. SPR is a classical and baseline physical layer model.

• Under the MPR model, up to M (M > 1) concurrently transmitted packets can be received successfully with non-zero probabilities as specified by a stochastic matrix Ξ defined as below:

Ξ                      ξ10 ξ11 ξ20 ξ21 ξ22 0 . . . . . . . . . . . . ξx 0 0 ξx 0 1 • • • • • • ξx 0 x 0 . . . . . . . . . . . . . . . ξM0 ξM1 • • • • • • ξMM 1 0 • • • • • • • • • 0                      (2.1)
where ξx 0 k 0 (k 0 ≤ x 0 ≤ M) is the probability of having k 0 successful packets among x 0 transmitted packets in a slot. Ξ is referred to as the reception matrix. The last two decades have witnessed an increasing prevalence of MPR technologies such as CDMA and MIMO. Mathematically, the SPR model can be regarded as a degenerated MPR model with M = 1 and

Ξ =         0 1 0 1 0 . . . . . . 1 0         .
The random access process operates as follows: FSA organises time-slots with each frame containing a number of consecutive time-slots. Each user is allowed to randomly and independently choose a time-slot to send his packet at most once per frame. More specifically, suppose the length of frame t is equal to L t , then in the beginning of frame t each user generates a random number R and selects the (R mod L t )-th time-slot in frame t to transmit his packet. Note that unsuccessful packets in the current frame are retransmitted in the next frame with the constant persistence probability p while newly generated packets are transmitted in the next frame following their arrivals with probability one.

For notation convenience, we use FSA-SPR and FSA-MPR to denote the FSA system operating on the SPR and MPR models, respectively.

Traffic model

Let random variable N t denote the total number of new arrivals during frame t and denote by A tl the number of new arrivals in time-slot l in frame t where l = 1, 2, • • • , L t . Assume that (A tl ) are independent and identically Poisson distributed random variables with probability distribution:

P{A tl = u} = Λ u (u ≥ 0) (2.2)
such that the expected number of arrivals per time-slot Λ = ∑ ∞ 1 uΛ u is finite. Then as N t = ∑ L t l=1 A tl , the distribution of N t , defined as {λ t (n)} n≥0 , also follows Poisson distribution with the expectation N t = L t Λ.

Packet success probability

The process of randomly and independently choosing a time-slot in a frame to transmit packets can be cast into a class of problems that are known as occupancy problems, or bins and balls problem [START_REF] Johnson | Urn models and their application: an approach to modern discrete probability theory[END_REF]. Specifically, consider the setting where a number of balls are randomly and independently placed into a number of bins, the classical occupance problem studies the maximum load of an individual bin.

In our context, time-slots and packets to be transmitted in a frame can be cast into bins and balls, respectively. Denote by Y t the random variable for the number of packets to be transmitted in frame t. Given Y t = ĥ in frame t and the frame length L t , the number x 0 of packets sent in one time-slot, referred as to occupancy number, is binomially distributed with parameters ĥ and 1 L t :

B ĥ, 1 L t (x 0 ) = ĥ x 0 ( 1 L t ) x 0 (1 - 1 L t ) ĥ-x 0 . (2.3)
Applying the distribution of equation (2.3) to all L t slots in the frame, we can get the expected value b(x 0 ) of the number of time-slots with occupance number x 0 in a frame as follows:

b(x 0 ) = L t B ĥ, 1 L t (x 0 ) = L t ĥ x 0 ( 1 L t ) x 0 (1 - 1 L t ) ĥ-x 0 . (2.4)
We further derive the probability that a packet is transmitted successfully under both SPR and MPR.

Packet success probability of FSA-SPR

In FSA-SPR, the number of successfully received packets equals that of time-slots with occupance number

x 0 = 1. Following the result of [START_REF] Vogt | Efficient object identification with passive RFID tags[END_REF], we can obtain the probability that under SPR there exist exactly k successful packets among ĥ transmitted packets in the frame, denoted by ξ SPR ĥk , as follows:

ξ SPR ĥk =                      ( L t k )( ĥ k )k!G(L t -k, ĥ-k) L ĥ t , 0 < k < min( ĥ, L t ) ( L t ĥ ) ĥ! L ĥ t , k = ĥ ≤ L t 0, k > min( ĥ, L t ) 0, k = L t < ĥ (2.5)
where

G(V, w) = V ŵ + ŵ ∑ t=1 (-1) t t-1 ∏ j=0 [( ŵ -j)(V -j)](V -t) ŵ-t 1 t!
with V L tk and ŵ ĥk.

Consequently, the expected number of successfully received packets in one frame in FSA-SPR, denoted as

r SPR h , is r SPR ĥ = min( ĥ,L t ) ∑ k=1 kξ SPR ĥk = b(1). (2.6)

Packet success probability of FSA-MPR

Let occupancy numbers x l and k l be the number of transmitted packets and successful packets in the lth time-slot, respectively, where l = 1, 2, • • • , L t . The probability that k packets are received successfully among ĥ transmitted packets in the frame, denoted by ξ MPR ĥk , can be expressed as

ξ MPR ĥk = ∑ ∑ l x l = ĥ ∑ ∑ l k l =k ∏ ξx l k l (2.7)
We can further derive the expected number of successfully received packets in one frame as

r MPR ĥ = ĥ ∑ k=1 kξ MPR ĥk = L t ĥ ∑ x 0 =1 x 0 ∑ k 0 =1 B ĥ, 1 L t (x 0 )k 0 ξx 0 k 0 . (2.8)
In the subsequent analysis, to make the presentation concise without introducing ambiguity, we use ξ ĥk to denote ξ SPR ĥk in FSA-SPR and ξ MPR ĥk in FSA-MPR. The notations used in the chapter are summarized in Table 2.1.

Main results

To streamline the presentation, we summarize the main results in this section and give the detailed proof and analysis in the subsequent sections that follow.

Aiming at studying the stability of FSA, we decompose our global objective into the following three questions, all of which are of fundamental importance both on the theoretical characterisation of FSA performance • Q2: When is the stability region maximised?

• Q3: How does FSA behave in the instability region?

Before answering the questions, we first introduce the formal definition of stability employed by Ghez et al.

in [START_REF] Ghez | Stability properties of slotted Aloha with multipacket reception capability[END_REF].

Define by random variable X t the number of backlogged packets in the system at the start of frame t. The discrete-time process (X t ) t≥0 can be seen as a homogeneous Markov chain.

Definition 2.1. An FSA system is stable if (X t ) t≥0 is ergodic and unstable otherwise.

By Definition 2.1, we can transform the study of stability of FSA into investigating the ergodicity of the backlog Markov chain. The rationality of this transformation is two-fold. One interpretation is the property of ergodicity that there exists a unique stationary distribution of a Markov chain if it is ergodic. The other can be interpreted from the nature of ergodicity that each state of the Markov chain can recur in finite time with probability 1.

From an engineering perspective, if FSA is stable, then the number of backlogs in the system will reduce overall; otherwise, it will increase as the system operates.

We then establish the following results characterizing the stability region and demonstrating the behavior of the Markov chain in nonergodicity regions under both SPR and MPR.

Results for FSA-SPR

Denote by i and ĥ the value of the number of backlogs and sent packets in frame t and α ĥ L t . Recall the definitions of X t and Y t , we can suppose that X t = i and Y t = ĥ. Theorem 2.1. Under FSA-SPR, consider an irreducible and aperiodic backlog Markov chain (X t ) t≥0 with nonnegative integers. When i→∞, we have 1 1. The system is always stable if Λ < αe -α and L t = Θ( ĥ). Specially, α = 1 maximizes the stability region2 and also the stable throughput.

2. The system is unstable under each of the following three conditions: (1) L t = o( ĥ); (2) L t = O( ĥ); (3) L t = Θ( ĥ) and Λ > αe -α .

Remark. Theorem 2.1 answers the first two questions and can be interpreted as follows:

• When L t = o( ĥ), i.e., the number of sent packets ĥ is far larger than the frame length L t , a packet experiences collision with high probability (w.h.p.), thus increasing the backlog size and destabilising the system;

• When L t = O( ĥ), i.e., the number of sent packets ĥ is far smaller than the frame length, a packet is transmitted successfully w.h.p.. However, the expected number of successful packets is still significantly less than that of new arrivals in the frame. The system is thus unstable.

• When L t = Θ( ĥ), i.e., ĥ has the same order of magnitude with the frame length, the system is stable when the backlog can be reduced gradually, i.e., when the expected arrival rate is less than the transmission success rate.

It is well known that an irreducible aperiodic Markov chain falls into one of three mutually exclusive classes: positive recurrent, null recurrent and transient. So, our next step after deriving the stability conditions is to show whether the backlog Markov chain in the instability region is transient or recurrent, which answers the third question. Theorem 2.2. With the same notations as in Theorem 2.1, (X t ) t≥0 is always transient in the instability region, i.e., under each of the following three conditions: (1) L t = o( ĥ); (2) L t = Θ( ĥ) and Λ > αe -α ; (3) L t = O( ĥ).

Remark. If a state of a Markov chain is transient, then the probability of returning to itself for the first time in a finite time is less than 1. Hence, Theorem 2.2 implies that once out of the stability region, the system is not guaranteed to return to stable state in finite time, that is, the number of backlogs will increase persistently. 1 For two variables X, Y, we use the following asymptotic notations:

• X=o(Y) * if 0 ≤ X Y ≤θ 0 , as Y→∞, where constant θ 0 ≥0; • X=o(Y) if X Y = 0, as Y → ∞; • X=O(Y) if X Y = ∞, as Y → ∞; • X=Θ(Y) if θ 1 ≤ X Y ≤θ 2 , as Y→∞, where constants θ 2 ≥ θ 1 > 0.

Results for FSA-MPR

Theorem 2.3. Under FSA-MPR, using the same notations as in Theorem 2.1, we have 1. The system is always stable if L t = Θ( ĥ) and Λ < ∑ M x 0 =1 e -α α x 0 x 0 ! ∑ x 0 k 0 =1 k 0 ξx 0 k 0 . Specially, let α * denote the value of α that maximises the upper bound of stability region, it holds that α * = Θ(M).

2.

The system is unstable under each of the following conditions: (1) L t = o( ĥ1-1 ) where 0 < 1 ≤ 1;

(2) L t = O( ĥ);

(3) Λ > α and L t = Θ( ĥ).

Remark.

Comparing the results of Theorem 2.3 to Theorem 2.1, we can quantify the performance gap between FSA-SPR and FSA-MPR in terms of stability. For example, when α = 1, the stability region is maximised in FSA-SPR with

Λ < e -1 , while the upper bound of the stability region in FSA-MPR is e -1 ∑ M

x 0 =1 1 (x 0 -1)! . Note that for M > 2, it holds that

1 + 1 + 1 2 < M ∑ x 0 =1 1 (x 0 -1)! < 1 + 1 + M ∑ x 0 =1 1 x 0 (x 0 + 1) < 2 + M ∑ x 0 =1 1 x 0 - 1 x 0 + 1 = 3 - 1 M + 1 .
The upper bound of the stability region of FSA-MPR when α = 1 is thus between 2.5 and 3 times the maximum stability region of FSA-SPR. And hence the maximum upper bound of the stability region of FSA-MPR achieved when

α * = Θ(M)
is far larger than that of FSA-SPR.

Theorem 2.4. With the same notations as in Theorem 2.3, (X t ) t≥0 is transient under each of the following three conditions:

(1) L t = o( ĥ1-1 ); (2) L t = O( ĥ); (3) Λ > α and L t = Θ( ĥ).

Remark. Theorem 2.4 demonstrates that despite the gain on the stability region size of FSA-MPR over FSA-SPR, their behaviors in the unstable region are essentially the same.

Stability Analysis of FSA-SPR

In this section, we will analyse the stability of FSA-SPR and prove Theorem 2.1 and 2.2.

Characterising backlog Markov chain

As mentioned in Sec. 2.4, we characterize the number of the backlogged packets in the system at the beginning of frame t as a homogeneous Markov chain (X t ) t≥0 . We assume that X t = i and Y t = ĥ. Denote by Z t the random variable for the number of retransmitted packets in frame t. Since the transmitted packets in frame t consists of the new arrivals during frame t -1 and the retransmitted packets in frame t, we have

Y t = Z t + N t-1 .
(2.9)

Suppose w new packets arrive in frame t -1 and h out of iw backlogs are retransmitted in frame t of which the probability is as follows:

B i-w (h) i -w h p h (1 -p) i-w-h .
As a consequence, the number of packets transmitted in frame t is ĥ = w + h.

We now calculate the one-step transition probability as a function of ξ ĥk , retransmission probability p and {λ t (n)} n≥0 . Denote by P is = P{X t+1 = s|X t = i} the one-step transition probability, we can derive the following results:

1) For i = 0:

P 00 = λ t (0), P 0s = λ t (s), s ≥ 1,
2) for i ≥ 1:

                       P i,i-s = i ∑ w=0 λ t-1 (w) i-w ∑ h={s-w} + B i-w (h) min(L, ĥ)-s ∑ n=0 λ t (n)ξ ĥ,n+s , 1 ≤ s ≤ i, P i,i = λ t (0) λ t-1 (0)B i (0) + i ∑ w=0 λ t-1 (w) i-w ∑ h=0 B i-w (h)ξ ĥ,0 + i ∑ w=0 λ t-1 (w) i-w ∑ h=0 B i-w (h) min( ĥ,L) ∑ n=1 λ t (n)ξ ĥn , P i,i+s = i ∑ w=0 λ t-1 (w) i-w ∑ h=0 B i-w (h) min( ĥ,L) ∑ n=0 λ t (n + s)ξ ĥn , s ≥ 1, (2.10)
where {s -w} + = max{sw, 0}.

The rationale for the calculation of the transition probability is explained as follows:

• When i = 0, i.e., there are no backlogs in the frame, the backlog size remains zero if no new packets arrive and increases by s if s new packets arrive in the frame.

• When i > 0, we have three possibilities, corresponding to the cases where the backlog size decreases, remains unchanged and increases, respectively:

-The state 1≤s≤min( ĥ, L t ) corresponds to the case where the backlog size decreases by s when n ≤ min( ĥ, L t )-s new packets arrive but n + s backlogged packets are received successfully.

-The backlog size remains unchanged if either of two following events happens: (a) no new packets are generated and either no backlogged packets are transmitted or all the transmitted backlogged packets fail; (b) n ≤ min( ĥ, L t ) new packets arrive but n backlogged packets are successfully received.

-The backlog size increases when the number of successful packets is less than that of new arrivals.

In order to establish the ergodicity of the backlog Markov chain (X t ) t≥0 , it is necessary to ensure (X t ) t≥0 is irreducible and aperiodic. To this end, we conclude this subsection by providing the sufficient conditions on {λ t (n)} for the irreducibility and the aperiodicity of (X t ) t≥0 as

0 < λ t (n) < 1, ∀ n ≥ 0. (2.11)
We would like to point out that most of traffic models can satisfy (2.11). Throughout the chapter, it is assumed that (2.11) holds and hence (X t ) t≥0 is irreducible and aperiodic.

Stability analysis

Recalling Definition 2.1, to study the stability of FSA, we need to analyse the ergodicity of the backlog Markov chain (X t ) t≥0 . We first define the drift and then introduce two auxiliary lemmas which will be useful in the ergodicity demonstration.

Definition 2.2. The drift D i of the backlog Markov chain (X t ) t≥0 at state X t = i where i ≥ 0 is defined as

D i = E[X t+1 -X t |X t = i].
(2.12)

Lemma 2.1 ([31]). Given an irreducible and aperiodic Markov chain (X t ) i≥0 having nonnegative integers as state space with the transition probability matrix [START_REF] Kaplan | A sufficient condition for nonergodicity of a Markov chain[END_REF]). Under the assumptions of Lemma 2.1, (X t ) t≥0 is not ergodic, if there exist some integer Q ≥ 0 and

P = {P is }, (X t ) t≥0 is ergodic if for some integer Q ≥ 0 and constant 0 > 0, it holds that 1. |D i | < ∞, for i ≤ Q, 2. D i < -0 , for i > Q. Lemma 2.2 ([
some constants B ≥ 0, c ∈ [0, 1] such that 1. D i > 0 for all i ≥ Q, 2. φ i -∑ s P is φ i ≥ -B(1 -φ) for all i ≥ Q, φ ∈ [c, 1].
Armed with Lemma 2.1 and Lemma 2.2, we start to prove Theorem 2.1.

Proof of Theorem 2.1. In the proof, we first explicitly formulate the drift defined by (2.12) and then study the ergodicity of Markov chain based on drift analysis.

Denote by random variable C t the number of successful transmissions in frame t, we have

X t+1 -X t = N t -C t .
Recall (2.12), it then follows that

D i = E[N t -C t |X t = i] = N t -E[C t |X t ]. (2.13)
Since all new arrivals and unsuccessful packets in frame t -1 are transmitted in frame t with probability one and p, respectively, we have

P{C t = k|X t = i, N t-1 = w, Z t = h} = ξ SPR ĥk ,
for 0 ≤ k ≤ min( ĥ, L). Recall (2.6), we have

E[C t |X t = i] = i-w ∑ h=0 B i-w (h)E[C t = k|X t = i, N t-1 = w] = i ∑ w=0 λ t-1 (w) i-w ∑ h=0 B i-w (h)r SPR ĥ .
(2.14) Following (2.13) and (2.14), we obtain the value of the drift as follows:

D i = N t - i ∑ w=0 λ t-1 (w) i-w ∑ h=0 B i-w (h)r SPR ĥ .
(2.15)

After formulating the drift, we then proceed by two steps.

Step 1: L t = Θ( ĥ) and Λ < αe -α .

In this step, we intend to corroborate that the conditions in Lemma 2.1 can be satisfied if L t = Θ( ĥ) and

Λ < αe -α . We first show that |D i | is finite. This is true for i ≤ Q since |D i | < max{N t , i ∑ w=0 λ t-1 (w) i-w ∑ h=0 B i-w (h)r SPR ĥ } < max{N t , min{L t , i ∑ w=0 λ t-1 (w) i-w ∑ h=0 (w + h)B i-w (h)}} < max{N t , min{L t , i ∑ w=0 wλ t-1 (w) + i ∑ w=0 (i -w)pλ t-1 (w)}} < max{N t , min{L t , (1 -p)λ + ip}}. (2.16)
Next, to derive the limit of D i , we start with the following lemma which is proved in Appendix 2.10.1.

Lemma 2.3. If r SPR ĥ has a limit r, then it holds that

lim i→∞ i ∑ w=0 λ t-1 (w) i-w ∑ h=0 B i-w (h)r SPR ĥ = r.
Following Lemma 2.3, we have where α ĥ L t . It thus holds that lim i→∞ D i < -0 with 0 = αe -α -Λ 2 since both α and Λ are constants when L t = Θ( ĥ) and Λ < αe -α .

lim i→∞ D i = N t -lim ĥ→∞ r SPR ĥ = lim ĥ→∞ L t Λ - ĥ 1 1 L t • 1 - 1 L t ĥ-1 = L t (Λ -αe -α ), (2.17
It then follows from Lemma 2.1 that (X t ) t≥0 is ergodic. Specially, when α = 1, the system stability region is maximized, i.e., Λ < e -1 .

Step 2: L t =o( ĥ) or L t =O( ĥ) or L t =Θ( ĥ) and Λ>αe -α .

In this step, we prove the instability of (X t ) t≥0 by applying Lemma 2.2. Taking into consideration the impact of different relation between L t and ĥ on the limit of D i . With (2.17), the following results hold for ĥ → ∞:

• Λlim α→∞ αe -α = Λ > 0, when L t = o( ĥ),

• Λlim α→0 αe -α = Λ > 0, when L t = O( ĥ),

• Λαe -α > 0, when L t = Θ( ĥ) and Λ > αe -α .

Consequently, we have lim i→∞ D i > 0, which proves the first condition in Lemma 2.2.

Next, we will validate the second condition of Lemma 2.2 in two cases according to the probable relationship between ĥ and i, i.e., ĥ = o(i) and ĥ = Θ(i).

Note that the second condition apparently holds for φ=0 and φ=1, we thus focus on the remaining value of φ, i.e., φ ∈ (c, 1). Moreover, given ĥ, P i,i-s in (2.10) can also be expressed as

P i,i-s = ĥ ∑ w=0 λ t-1 (w)B i-w ( ĥ -w) ĥ-s ∑ n=0 λ t (n)ξ ĥ,n+s . (2.18)
Now, we start the proof with the above arms.

Case 1: ĥ = o(i).

Given ĥ = o(i), we can derive the result as follows:

∞ ∑ s=0 φ s P is = i-ĥ-1 ∑ s=0 φ s P is + i ∑ s=i- ĥ φ s P is + ∞ ∑ s=i+1 φ s P is ≤ φ i+1 + i ∑ s=i- ĥ φ s ĥ ∑ w=0 λ t-1 (w)B i-w ( ĥ -w) • ĥ+s-i ∑ n=0 λ t (n)ξ ĥ,n+i-s ≤ φ i+1 + i ∑ s=i- ĥ φ s ĥ ∑ w=0 B i-w ( ĥ -w) ≤ φ i+1 + ĥe -ip 2 (1-ĥ ip ) 2 φ i-ĥ ≤ φ i , as i → ∞, (2.19)
where we use the Chernoff's inequality to bound the cumulative probability of B i-w ( ĥw). Therefore, the second condition of Lemma 2.2 holds when ĥ = o(i).

Case 2: ĥ = Θ(i).

In this case, we need to distinguish the three instability regions. Without loss of generality, we assume that ĥ = βi where constant β ∈ (0, 1].

(1) L t = o( ĥ).

When L t = o( ĥ), it also holds that L t = o(i) and that at most L t -1 packets are successfully received, we thus have

∞ ∑ s=0 φ s P is = i-L t -1 ∑ s=0 φ s P is + i ∑ s=i-L t φ s P is + ∞ ∑ s=i+1 φ s P is ≤ φ i+1 + L t φ i-L t ≤ φ i + B(1 -φ), as i → ∞, (2.20)
for any positive constant B.

(2) L t = Θ( ĥ).

The key steps we need in this case are to obtain upper bounds of ξ SPR ĥk and the arrival rate in a new way. To this end, we first recomputed ξ SPR ĥk when L t = Θ( ĥ) as follows:

                     ξ SPR ĥk = ( L k )k!(L t -k) ĥ-k L ĥ t - ( L t k+1 )(k+1)!(L t -k-1) ĥ-k-1 L ĥ t ≤ ( L t k )k! L k t (1 -k L t ) ĥ-k -(1 -k+1 L t ) ĥ-k ≤ ( L t k )k! L k t ≤ (1 -k/2 L t ) k/2 , k < min( ĥ, L t ),
ξ SPR ĥ ĥ = ( L t ĥ ) ĥ! L ĥ t ≤ (1 -ĥ/2 L t ) ĥ/2 ≤ (1α/2) ĥ/2 , k = ĥ ≤ L t .

(2.21)

The rationale behind the above inequalities is as follows: Given ĥ transmitted packets, the probability of exactly k successful packets equals the absolute value of the difference between the probability of at least k successful packets and that of at least k + 1 successful packets.

Next, we introduce an auxiliary lemma to bound the probability distribution of the arrival rate. When the number of new arrivals per slot A tl is Poisson distributed with the mean Λ, the number of new arrivals per frame N t (A tl and N t is formally defined in Sec. 2.3.) is also a Poisson random variable with the mean

N t = L t Λ > ĥe -α .
Lemma 2.4 ([33]). Given a Poisson distributed variable X with the mean µ, it holds that

Pr[X ≤ x] ≤ e -µ (eµ) x x x , ∀ x < µ, (2.22) Pr[X ≥ x] ≤ e -µ (eµ) x x x , ∀ x > µ.
(2.23)

In the case that L t = Θ(i), it holds that

N t = L t Λ > L 2 3
t , for the constant Λ and a large i. Consequently, applying (2.22) in Lemma 2.4, we have

P{N t ≤ L 2/3 t } ≤ e -λ (eλ) L 2/3 t (L 2/3 t ) L 2/3 t ≤ e -L 2/3 t ( LΛ L 2/3 t -1) L t Λ L 2/3 t L 2/3 t ≤ e L t Λ L 2/3 t -1 L 2/3 L t Λ -L 2/3 t ≤ 1 a L 2/3 t 1 , (2.24)
where a 1 e ΛL 1/3 t -1 ΛL 1/3 t 1, following the fact that e x > 1 + x, for ∀ x > 0.

Armed with (2.21) and (2.24) and noticing the fact that at most ĥ packets are successfully received, we start developing the proof and obtain the results as follows:

∞ ∑ s=0 φ s P is = i-ĥ-1 ∑ s=0 φ s P is + i ∑ s=i- ĥ φ s P is + ∞ ∑ s=i+1 φ s P is ≤ φ i+1 + i ∑ s=i- ĥ ĥ ∑ w=0 λ t-1 (w)B i-w ( ĥ -w) • ĥ+s-i ∑ n=0 λ t (n)(1 - n + i -s 2L t ) (n+i-s)/2 ≤ φ i+1 + i ∑ s=i- ĥ φ s ĥ ∑ n=0 λ t (n)(1 - n 2L t ) n 2 ≤ φ i+1 + i ∑ s=i- ĥ φ s L 2/3 t ∑ n=0 λ t (n) + ĥ ∑ n=L 2/3 t λ t (n)(1 - n 2L t ) n 2 ≤ φ i+1 + i ∑ s=i- ĥ φ s 1 a L 2/3 t 1 + (1 - 1 2L 1/3 t ) L 2/3 2 ≤ ( ĥ + 1) 1 a L 2/3 t 1 + e -L 1/3 t 4 φ i-ĥ + φ i+1 ≤ φ i + B(1 -φ), as i → ∞, (2.25)
for any positive constant B, where the last inequality holds for ( ĥ + 1) 1

a L 2/3 t 1 + e -L 1/3 t 4 ∼ Θ(ie -i 1/3 )→ 0 as i → ∞, while B(1 -φ) is positive constant.
Consequently, the second condition in Lemma 2.2 holds for Case 2. Next, we proceed with the proof for the third case.

(3) L t = O( ĥ).

When L t = O( ĥ), it also holds that L t = O(i) such that the expected number of new arrivals per frame

N t = L t Λ i.
φ s P is = i ∑ s=0 φ s P is + ∞ ∑ s=i+1 φ s P is ≤ i ∑ s=0 φ s s ∑ n=0 λ t (n) + φ i+1 ≤ i ∑ s=0 s ∑ n=0 λ t (n) + φ i+1 ≤ i + 1 (φ L t Λ i ) i φ i + φ i+1 ≤ φ i , as i → ∞, (2.27)
since φ is constant while L t Λ i → ∞ as i → ∞. Combining the analysis above, it follows Lemma 2.2 that the backlog Markov chain (X t ) t≥0 is unstable when L t = o( ĥ) or L t = O( ĥ) or L t = Θ( ĥ) and Λ > αe -α . And the proof of Algorithm 2.1 is thus completed.

System behavior in instability region

It follows from Theorem 2.1 that the system is unstable in the following three conditions: L t = o( ĥ); L t = O( ĥ); and L t = Θ( ĥ) but Λ > αe -α . Lemma 2.2, however, is not sufficient to ensure the transience of a Markov chain, we thus in this section further investigate the system behavior in the instability region, i.e., when (X t ) t≥0 is nonergodic. The key results are given in Theorem 2.2.

Before proving Theorem 2.2, we first introduce the following lemma [START_REF] Mertens | Necessary and sufficient conditions for recurrence and transience of Markov chains, in terms of inequalities[END_REF] on the conditions for the transience of a Markov chain. Lemma 2.5 ([34]). Let (X t ) t≥0 be an irreducible and aperiodic Markov chain with the nonnegative integers as its state space and one-step transition probability matrix P = {P is }. (X t ) t≥0 is transient if and only if there exists a sequence

{y i } i≥0 such that 1. y i (i ≥ 0) is bounded, 2. for some i ≥ N, y i < y 0 , y 1 , • • • , y N-1 , 3. for some integer N > 0, ∑ ∞ s=0 y s P is ≤ y i , ∀ i ≥ N.
Armed with Lemma 2.5, we now prove Theorem 2.2.

Proof of Theorem 2.2. The key to prove Theorem 2.2 is to show the existence of a sequence satisfying the properties listed in Lemma 2.5, so we first construct the following sequence (2.28) and then prove that it satisfies the required conditions.

y i = 1 (i + 1) θ , θ ∈ (0, 1).
(2.28)

It can be easily checked that {y i } satisfies the first two properties in Lemma 2.5.

Noticing that the sequence {φ i } in Lemma 2.2 satisfies the first two properties in Lemma 2.5 for 0 < φ < 1, and recall (2. [START_REF] Szpankowski | Stability conditions for some distributed systems: Buffered random access systems[END_REF]) and (2.27), we can conclude that (X t ) t≥0 is transient if ĥ = o(i) or L t = O(i). Therefore, we next proceed with ĥ = Θ(i) by distinguish two cases.

Case 1: L t = o( ĥ).

When ĥ = Θ(i), it also holds that L t = o(i). To streamline the complicated analysis in this case, we partition the region L t =o(i) into two parts, i.e., 1) L t =o((ln i) 4 ) * , and 2) L t =o(i) except part 1), i.e., the region [O((ln i) 4 ), o(i)].

• Part 1): L t =o((ln i) 4 ) * .

The result in this part is shown in the following lemma for the third property in Lemma 2.5. The proof is detailed in Appendix 2.10.2.

Lemma 2.6. If L t =o((ln i) 4 ) * , then (X t ) t≥0 is always transient.

• Part 2): L t =o(i) except part 1).

In this case, since a 1 > ln i and y i -

y i+1 = 1 (i+1) θ (1 -(1 -1 i+2 ) θ ) ≥ θ (i+1) θ (i+2)
where we use the fact that (1 -1 i+2 ) θ ≤ 1 -θ i+2 following Taylor's theorem, using (2.24) and (2.28) yields

∞ ∑ s=0 y s P is = i-L t ∑ s=0 y s P is + i ∑ s=i-L t +1 y s P is + ∞ ∑ s=i+1 y s P is ≤y i+1 + i ∑ s=i-L t +1 y s s ∑ n=0 λ t (n)(1 - n 2L t ) n 2 ≤ i ∑ s=i-L t +1 y s L 2 3 t ∑ n=0 λ n + s ∑ n=L 2 3 t +1 λ t (n)(1 - n 2L t ) n 2 + y i+1 ≤ i ∑ s=i-L t +1 y s 1 a L 2/3 t 1 + (1 - 1 2L 1/3 t ) L 2/3 t 2 + y i+1 ≤ L (i -L t + 2) θ 1 a L 2/3 t 1 + e -L 1/3 t 4 + 1 (i + 2) θ ≤ L t (i -L t + 2) θ (ln i) -(ln i) 8/3 + i -(ln i) 1/3 4 + 1 (i + 2) θ ≤i -4 + 1 (i + 2) θ ≤ 1 (i + 1) θ , as i → ∞.
(2.29)

Case 2: L t = Θ( ĥ).

In this case, the method to prove is similar with that used in (2.25). Recall (2.25) , we have

∞ ∑ s=0 y s P is = i-ĥ-1 ∑ s=0 y s P is + i ∑ s=i- ĥ y s P is + ∞ ∑ s=i+1 y s sP is ≤ y i+1 + i ∑ s=i- ĥ y s 1 a L 2/3 1 + (1 - 1 2L 1/3 ) L 2/3 2 ≤ ĥ + 1 (i -ĥ + 1) θ 1 a L 2/3 1 + e -L 1/3 4 + 1 (i + 2) θ ≤ 1 (i + 1) θ , as i → ∞.
(2.30)

Consequently, it follows Lemma 2.5 that the backlog Markov chain (X t ) t≥0 is transient in the instability region, which completes the proof of Theorem 2.2.

Stability Analysis of FSA-MPR

In this section, we study stability properties of FSA-MPR. Following a similar procedure as the analysis of FSA-SPR, we first establish conditions for the stability of FSA-MPR and further analyse the system behavior in the instability region.

Stability analysis

We employ Lemma 2.1 and Lemma 2.2 as mathematical base to study the stability properties of FSA-MPR, more specifically, in the proof of Theorem 2.3.

Proof of Theorem 2.3. We develop our proof in 3 steps.

Step 1: stability conditions.

In step 1, we prove the conditions for the stability of (X t ) t≥0 , i.e., Λ< ∑ M x 0 =1 e -α α x 0 x 0 ! ∑ x 0 k 0 =1 k 0 ξx 0 k 0 and L t = Θ( ĥ).

Similar to (2.15), the drift at state i of (X t ) t≥0 in FSA-MPR can be written as:

D i = N t - i ∑ w=0 λ t-1 (w) i-w ∑ h=0 B i-w (h)r MPR ĥ .
(2.31)

According to (2.16), D i is finite as shown in the following inequality:

|D i | < max{N t , (1 -p)N t + ip},
which demonstrates the first conditions in Lemma 2.1 for the ergodicity of (X t ) t≥0 .

Recall (2.8) and Lemma 2.3, we have

lim i→∞ D i = lim i→∞ N t - i ∑ w=0 λ t-1 (w) i-w ∑ h=0 B i-w (h)r MPR ĥ = lim ĥ→∞ L t Λ -ĥ ∑ x 0 =1 B ĥ, 1 L (x 0 ) x 0 ∑ k 0 =1 k 0 ξx 0 k 0 = L t Λ - M ∑ x 0 =1 e -α α x 0 x 0 ! x 0 ∑ k 0 =1 k 0 ξx 0 k 0 . (2.32)
Therefore, it holds that lim i→∞ D i < -0 if L t = Θ( ĥ) and

Λ < ∑ M x 0 =1 e -α α x 0 x 0 ! ∑ x 0 k 0 =1 k 0 ξx 0 k 0 R1 . It then follows from Lemma 2.1 that (X t ) t≥0 is ergodic with 0 = R1 -Λ 2 .
Step 2: α * = Θ(M).

In Step 2, we show that α * = Θ(M). Since the proof consists mainly of algebraic operations of function optimization, we state the following lemma proving Step 2 and detail its proof in Appendix 2.10.3. Lemma 2.7. Let α * denote the value of α that maximises the upper bound of the stability region, it holds that α * = Θ(M).

Step 3: instability region.

In Step 3, we prove the instability region of (X i ) i≥0 by applying Lemma 2.2.

When L t = o( ĥ), recall (2.32), we have ĥ

∑ x 0 =1 B ĥ, 1 L (x 0 ) x 0 ∑ k 0 =1 k 0 ξx 0 k 0 = M ∑ x 0 =1 B ĥ, 1 L (x 0 ) x 0 ∑ k 0 =1 k 0 ξx 0 k 0 ≤ M ∑ x 0 =1
x 0 B ĥ, 1 L (x 0 ) → 0, as ĥ → ∞, since lim ĥ→∞ B ĥ, 1 L (x 0 ) = 0 for a finite M.

Moreover, for L t = O( ĥ), it can be derived from (2.32)

that ∑ M x 0 =1 e -α α x 0 x 0 ! ∑ x 0 k 0 =1 k 0 ξx 0 k 0 →0 since α → 0 as ĥ → ∞.
Furthermore, according to the analysis in the first step, we know that lim i→∞ D i > 0, if the conditions in the first step are not satisfied.

Additionally, in the analysis of FSA-SPR system, we have proven that if ĥ = o(i) or L t =o( ĥ) or L t =O( ĥ), the Markov chain (X t ) t≥0 is always unstable, independent of ξ ĥk . Noticing that ξ ĥk is the only difference between FSA-SPR and FSA-MPR, it thus also holds that (X t ) t≥0 is unstable under FSA-MPR in the three cases.

We next study the instability of FSA-MPR when L t = Θ( ĥ) and Λ > α. In this case, it holds that N t = L t Λ > ĥ such that

P{N t ≤ ĥ} ≤ 1 a ĥ 3 , (2.33)
where a 3 α Λ e Λ α -1 > 1. Note that the one-step transition probability P is in FSA-MPR can be obtained by replacing min( ĥ, L t ) with ĥ in (2.10).

Hence, recall (2.25), we have

∞ ∑ s=0 φ s P is = i ∑ s=0 φ s P is + ∞ ∑ s=i φ s P is ≤ i ∑ s=0 φ s ĥ ∑ n=0 λ t (n)ξ ĥ,n+i-s + φ i+1 ≤ 1 a βi 3 + φ i+1 ≤ φ i + B(1 -φ), as i → ∞, (2.34)
which proves the instability of FSA-MPR following Lemma 2.2 and also completes the proof of Theorem 2.3.

System behavior in instability region

It follows from Theorem 2.3 that the system is unstable under the following three conditions: L t = o( ĥ);

L t = O( ĥ); L t = Θ( ĥ) and Λ > α. In this subsection, we further investigate the system behavior in the instability region, i.e., when (X t ) t≥0 is nonergodic. The key results are given in Theorem 2.4, whose proof is detailed as follows.

Proof of Theorem 2.4. In the proof of Theorem 2.2, we have proven that when L t =O( ĥ) or ĥ=o(i), the Markov chain (X t ) t≥0 is always transient, we thus develop the proof for ĥ=Θ(i) by distinguishing two cases.

Case 1:

L t = o( ĥ1-1 ) * with 1 ∈ (0, 1].
In this case, it holds that L t = o(i 1-1 ) * for ĥ=Θ(i). As counterparts in FSA-SPR, we also partition the region into two parts, i.e., 1) L t =o((ln i) 4 ) * , and 2) L t =o(i) except part 1), i.e., the region [O((ln i) 4 ), o(i)].

Recall the proof of Lemma 2.6, it has been shown that (X t ) t≥0 is always transient, independent of ξ ĥk , meaning (X t ) t≥0 is also transient in FSA-MPR when L t =o((ln i) 4 ) * .

As a consequence, it is sufficient to show the transience of (X t ) t≥0 in part 2). The key step here is to obtain the upper bound of ξ ĥk . To this end, we first introduce the following auxiliary lemma.

Lemma 2.8 ([35]). Given ĥ packets, each packet is sent in a slot picked randomly among L t time-slots in frame t. If

ρ j = L t
e -ĥ/L t j! ( ĥ L t ) j remains bounded for ĥ, L t → ∞, then the probability P(m j ) of finding exactly m j time-slots with j packets can be approximated by the following Poisson distribution with the parameter ρ j ,

P(m j ) = e -ρ j ρ j m j m j ! . (2.35)
We next show that Lemma 2.8 is applicable to FSA-MPR when L t = o( ĥ1-1 ) * for a large enough ĥ. To that end, we verify the boundedness of ρ j , which is derived as

0 ≤ ρ j ≤ ĥj j!L j-1 t e ĥ ≤ ĥj j!L j-1 t • ( 1 j)!
( ĥ )

1 j ≤ ( 1 j)! j!L j-1 t , (2.36)
meaning that ρ j is bounded if j is finite.

Apparently, when L t = o( ĥ1-1 ) * , the probability of finding exactly m j time-slots with j packets in FSA-MPR can be approximated by the Poisson distribution with the parameter ρ j , following from Lemma 2.8 with

j = 1, 2, • • • , M.
Consequently, we can derive the probability ξ MPR 1→M that there are no slots with 1 ≤ j ≤ M packets as follows:

ξ MPR 1→M = e -(ρ 1 +ρ 2 +•••+ρ M ) .
(2.37)

Furthermore, since the event that all ĥ packets fail to be received has two probabilities, i.e., 1) there are no slots with 1 ≤ j ≤ M packets in the whole frame, and 2) there exists slots with 1 ≤ j ≤ M packets, but all of these packets are unsuccessful. As a result, it holds that ξ MPR ĥ0 ≥ ξ MPR 1→M . We thus can get the following inequalities:

ξ MPR ĥk ≤ 1 -ξ MPR ĥ0 ≤ 1 -e -(ρ 1 +ρ 2 +•••+ρ M ) ≤ 1 -e -Mρ M , k ≥ 1, (2.38)
where we use the fact that the probability of exact k ≥ 1 successfully received packets among ĥ packets is less than that of at least one packet received successfully in the first inequality. And the third inequality above follows from the monotonicity of ρ j when L = o( ĥ1-1 ) * , i.e.,

ρ M > ρ M-1 > • • • > ρ 2 > ρ 1 .
In addition, we can also derive the following results:

0 ≤ lim ĥ→∞ ĥ4 (1 -e -Mρ M ) ≤ lim ĥ→∞ e Mρ M -1 (1/ ĥ4 ) ≤ lim ĥ→∞ 1 ĥM 1 +5 4M!e ĥ 1 ≤ lim ĥ→∞ ∏ M-1+ 5 1 x=0 (M + 5 1 -x) 4M!e ĥ 1 ≤ 0, which means 1 -e -Mρ M ≤ 1 ĥ4 .
Using this inequality and recall (2.39), we have

∞ ∑ s=0 y s P is = i-L ∑ s=0 y s P is + i ∑ s=i-L+1 y s P is + ∞ ∑ s=i+1 y s P is ≤ y i+1 + i ∑ s=i-L+1 y s s ∑ n=0 λ t (n)ξ ĥ,n+i-s ≤ i ∑ s=i-L+1 y s L 2 3 ∑ n=0 λ n + s ∑ n=L 2 3 +1 λ t (n)ξ ĥ,n+i-s + y i+1 ≤ i ∑ s=i-L+1 y s 1 a L 2/3 1 + 1 ĥ4 + y i+1 ≤ L (i -L + 2) θ 1 a L 2/3 1 + 1 ĥ4 + 1 (i + 2) θ ≤ 2(βi) -3 + 1 (i + 2) θ ≤ 1 (i + 1) θ , as i → ∞.
(2.39)

Thus, according to Lemma 2.6, the backlog Markov chain (X t ) t≥0 is transient when

L t = o( ĥ1-1 ) * .
Case 2: L t = Θ( ĥ) and Λ > α.

In this case, we have N t = L t Λ > ĥ. Using similar reasoning as (2.34), we have

∞ ∑ s=0 y s P is ≤ βi + 1 a βi 3 + 1 (i + 2) θ ≤ 1 (i + 1) θ , as i → ∞.
Therefore, (X t ) t≥0 is also transient in this case and the proof of Theorem 2.4 is completed.

Discussion

In the sections above, we prove that the stability of FSA relies on the relationship between the frame size and the number of packets to be transmitted in this frame. In order to stabilize FSA systems, all the users should know the value of transmitted packets in the current frame. The state information, However, is imperfect in some scenarios such that the users do not have access to the value of packets to be transmitted in the frame.

Fortunately, we can get its approximate value.

Recall (2.9), because N t follows the Poisson distribution and Z t follows the binomial distribution which can be approximated as the Poisson distribution, Y t can also be approximated as a Poisson distributed random variable. According to Lemma 2.4, the value of Y t sharply concentrates around its expectation, we thus use the following E[Y t ] to approximate ĥ:

E[Y t ] = ip + (1 -p)E[N t-1 |X t = i] = ip + (1 -p) i ∑ w=0 we -L t-1 Λ (L t-1 Λ) w w! . (2.40)
As a result, we can set the frame size following the control algorithm as follows:

L t = c 1 ip + (1 -p) i ∑ w=0 we -L t-1 Λ (L t-1 Λ) w w! , ( 2.41)
L 0 = c 1 ϑ, (2.42)
where X 0 =ϑ means the initial number of packets in the system, and c 1 = 1 for FSA-SPR and c 1 = 1 α * for FSA-MPR.

By the above control algorithm, the frame size L t only depends on the value of backlog population size i, so the original problem is translated to estimate the number of backlogs X t , i.e., i. Fortunately, there exist several estimation approaches which exploit the channel feedback, such as the probability of a idle or collision slot, and the number of idle or collision slots.

According to the requirement on the estimation accuracy, we can select a rough estimator or an accurate estimator. Since ip ≤ L t ≤ i in (2.41), we can estimate X t roughly so that the estimate Xt = Θ(X t ) in very short time, more specifically, in log(X t ) or log log(X t ) slots [START_REF] Chen | Understanding RFID counting protocols[END_REF]. While if the accurate result is required, we can use the additive estimator as in [START_REF] Wu | Fasa: Accelerated S-ALOHA using access history for event-driven M2M communications[END_REF] and Kalman filter-based estimator as in our another work [START_REF] Yu | From static to dynamic tag population estimation: An extended Kalman Filter perspective[END_REF] to estimate the value of X t and update L t .

Numerical Results

In this section, we conduct simulations via MATLAB to verify our theoretical results by illustrating the evolution of the number of backlogs in each frame under different parameters with the following default settings: the initial number of backlogs X 0 =10 4 , the simulation duration t max =100, o( ĥ)≤0.01 ĥ, O( ĥ) ≥ 100 ĥ, α = 1 in FSA-SPR and α = α * ∈ ( M-1 e , M) in FSA-MPR when L t = Θ( ĥ). To simulate FSA, each user first generates a random number among [0, L t -1] uniformly and responds in the corresponding slot. And all results are obtained by taking the average of 100 trials.

Stability properties of FSA

FSA-SPR systems. We start by investigating numerically the stability properties of FSA-SPR. As stated in Theorem 2.1, if Λ < 1 e , the system is stable and unstable otherwise when α = 1, we thus set the expected arrival rate per slot Λ to 0.3 and 0.37 for the analysis of stability and instability for L t = ĥ, respectively. Moreover, we set a small Λ = 0.01 to analyze the instability for the cases L t = o( ĥ) and L t = O( ĥ).

As shown in Fig. 2.1(a), for the case L t = ĥ, the number of backlogs decreases to zero at a rate in proportion to the retransmission probability if Λ < 1 e , while increasing gradually otherwise. This is due to the nature of FSA that frame size varies with the number of the sent packets to maximize the throughput per slot. Moreover, Fig. 2.1(b) and Fig. 2.1(c) illustrate the instability when L t = o( ĥ) and L t = O( ĥ). The numerical results is in accordance with the analytical results on FSA-SPR in Theorem 2.1.

FSA-MPR System.

We then move to the FSA-MPR exploiting MPR model as in [START_REF] Zeng | Neighbor discovery in wireless networks with multipacket reception[END_REF] with M = 10. Recall Lemma 2.7, it can be derived that α * =10/1.37 and the maximum stability region is R1 = 5.814. We thus set Λ = 5 and Λ = 5.9 for the stability analysis in the case L t = ĥ/α * . From Fig. 2.2, we can see that the numerical results is in accordance with the analytical results on FSA-MPR.

Comparison under different frame sizes

Here we evaluate the performance difference when the frame size deviates from its optimum value that α = ĥ L t = 1 in FSA-SPR and α = α * in FSA-MPR. To that end, we set Λ = 0.3 for FSA-SPR and Λ = 5 for FSA-MPR. As shown in Fig. 2.3 and Fig. 2.4, the performance degrades significantly when the frame size is not optimal because the throughput is reduced in this case.

Comparison between FSA-SPR and FSA-MPR

We further compare the performance of FSA-SPR and FSA-MPR. To that end, for both FSA-SPR and FSA-MPR, we set Λ = 0.3 and L t = ĥ where ĥ is the number of backlogs in FSA-SPR maximizing the throughput of FSA-SPR. Although this setting is not optimal to FSA-MPR, we can also see from Fig. 2.5 that FSA-MPR remarkably outperforms FSA-SPR.

Conclusion

In this chapter, we have studied the stability of FSA-SPR and FSA-MPR by modeling the system backlog as a Markov chain. By employing drift analysis, we have obtained the closed-form conditions for the stability of FSA and shown that the stability region is maximised when the frame length equals the number of sent packets in FSA-SPR and the upper bound of stability region is maximised when the ratio of the number of sent packets to frame length equals in an order of magnitude the maximum multipacket reception capacity in FSA-MPR. Furthermore, to characterise system behavior in the instable region, we have mathematically demonstrated the existence of transience of the Markov chain. In addition, we conduct the numerical analysis to verify the theoretical results. Our results provide theoretical guidelines on the design of stable FSA-based protocols in practical applications such as RFID systems and M2M networks.

Proofs

Proof of Lemma 2.3

Proof. We take into account two cases depending on whether r is finite.

Case 1: r = +∞.
Recall the definition of the limit, fix M > 0 and pick a w * ≥0 such that λ t-1 (w * ) =0, there exists an integer Q such that r SPR ĥ >M for all ĥ≥Q. Then fix such Q, we have

i ∑ w=0 λ t-1 (w) i-w ∑ h=0 B (i-w) (h)r SPR h+w > λ t-1 (w * ) i-w * ∑ h=0 B (i-w * ) (h)r SPR h+w * > λ t-1 (w * ) • i-w * ∑ h=Q B (i-w * ) (h)r SPR h+w * > Mλ t-1 (w * ) i-w * ∑ h=Q B (i-w * ) (h)
for all i ≥ Q + w * , since for any fixed Q and w * it holds that

lim i→∞ i-w * ∑ h=Q B (i-w * ) (h) = 1. (2.43)
Case 2: r < +∞.

For i > 2Q, after some algebraic operations, we get

i ∑ w=0 λ t-1 (w) i-w ∑ h=0 B (i-w) (h)r SPR h+w = Q ∑ w=0 λ t-1 (w) i-w ∑ h=0 B (i-w) (h)r SPR h+w + i ∑ w=Q+1 λ t-1 (w) i-w ∑ h=0 B (i-w) (h)r SPR h+w = Q ∑ w=0 λ t-1 (w) 2Q-w ∑ h=0 B (i-w) (h)r SPR h+w + Q ∑ w=0 λ t-1 (w) i-w ∑ h=2Q-w+1 B (i-w) (h)r SPR h+w + i ∑ w=Q+1 λ t-1 (w) i-w ∑ h=0 B (i-w) (h)r SPR h+w . (2.44)
Now we proceed by calculate three items in the right hand side of the last equality in (2.44).

• In the the first item, since 0 ≤ r h+w ≤ h + w ≤ 2Q and 0 < λ t-1 (w) < 1 and (2.43) holds, we have

∑ Q w=0 λ t-1 (w) ∑ 2Q-w h=0 B (i-w) (h)r SPR h+w → 0, as i → ∞. • In the second item, since it holds that h + w > Q and (2.43), we have lim i→∞ Q ∑ w=0 λ t-1 (w) i-w ∑ h=2Q-w+1 B (i-w) (h)r SPR h+w = lim i→∞ r Q ∑ w=0 λ t-1 (w) i-w ∑ h=2Q-w+1 B (i-w) (h) = r Q ∑ w=0 λ t-1 (w).
• In the third item, since it holds that h + w > Q, we have

lim i→∞ i ∑ w=Q+1 λ t-1 (w) i-w ∑ h=0 B (i-w) (h)r SPR h+w = lim i→∞ r i ∑ w=Q+1 λ t-1 (w) i-w ∑ h=0 B (i-w) (h) = r i ∑ w=Q+1 λ t-1 (w).
Consequently, we can get

lim i→∞ i ∑ w=0 λ t-1 (w) i-w ∑ h=0 B (i-w) (h)r SPR h+w = r Q ∑ w=0 λ t-1 (w) + r i ∑ w=Q+1 λ t-1 (w) = r, (2.45)
which completes the proof of Case 2 and the lemma as well.

Proof of Lemma 2.6

Proof. Recall the definition of transition probability, we can get the following equivalent:

∑ s y s P is ≤ y i ⇐⇒ min(i,L t ) ∑ s=1 (y i-s -y i )P i,i-s + ∞ ∑ s=1 (y i+s -y i )P i,i+s ≤ 0.
With the following definitions:

             f (i) = (i + 1) θ min(i,L t ) ∑ s=1 1 (i + 1 -s) θ - 1 (i + 1) θ i ∑ w=0 λ t-1 (w) i-w ∑ h={s-w} + B i-w (h) min(L t , ĥ)-s ∑ n=0 λ t (n)ξ ĥ,n+s , g (i) = (i + 1) θ ∞ ∑ s=1 1 (i + 1 + s) θ - 1 (i + 1) θ i ∑ w=0 λ t-1 (w) i-w ∑ h=0 B i-w (h) min( ĥ,L t ) ∑ n=0 λ t (n + s)ξ ĥn , (2.46) we have ∑ s y s P is ≤ y i ⇐⇒ f (i) + g (i) ≤ 0.
Moreover, the drift can be rewritten as

D i = - min(i,L t ) ∑ s=1 sP i,i-s + ∞ ∑ s=1 sP i,i+s = f (i) + g(i),
where f (h) and g(h) are defined as

             f (i) = min{i,L t } ∑ s=1 s i ∑ w=0 λ t-1 (w) i-w ∑ h={s-w} + B i-w (h) min(L t , ĥ)-s ∑ n=0 λ t (n)ξ ĥ,n+s , g(i) = ∞ ∑ s=1 s i ∑ w=0 λ t-1 (w) i-w ∑ h=0 B i-w (h) min( ĥ,L t ) ∑ n=0 λ t (n + s)ξ ĥn . (2.47) Therefore, (X t ) t≥0 is transient if it holds that lim i→∞ [ f (i) + g (i) + θD i ] = 0.
(2.48)

Noticing that D i = f (i) + g(i), we prove (2.48) by showing that: (1) lim i→∞ [ f (i) + θ f (i)] = 0; (2) lim i→∞ [g (i) + θg(i)] = 0. We first prove lim i→∞ [ f (i) + θ f (i)] = 0.
From (2.46) and (2.47), we get

f (i) + θ f (i) = (i + 1) min(i,L t ) ∑ s=1 ( i + 1 i + 1 -s ) θ -1 - θs i + 1 • i ∑ w=0 λ t-1 (w) i-w ∑ h={s-w} + B i-w (h) L t -s ∑ n=0 λ t (n)ξ ĥ,n+s , (2.49)
which is nonnegative since

( i + 1 i + 1 -s ) θ -1 - θs i + 1 > 0, ∀1 ≤ s ≤ i.
and thus

f (i) + θ f (i) ≤ (i + 1) L t ∑ s=1 ( i + 1 i + 1 -s ) θ -1 - θs i + 1 . Given 0 < v ≤ L t < i, define m i (v) as m i (v) = i + 1 v 2 ( i + 1 i + 1 -v ) θ -1 - θ v
which is positive and nondecreasing in v for i ≥ 1, we get

m i i + 1 2 ≤ 1 i + 1 [4(2 θ -1) -2θ] A i + 1 ,
where A is a positive constant only depending on θ. Thus,

f (i) + θ f (i) ≤ L t ∑ s=1 s 2 m i (s) ≤ A i + 1 L t ∑ s=1 s 2 ≤ AL t (L t + 1)(2L t + 1) 6(i + 1) . Since lim i→∞ AL t (L t +1)(2L t +1) 6(i+1) = 0 when L t = o((ln i) 4 ) * , it thus holds that lim i→∞ [ f (i) + θ f (i)] = 0.
We then prove lim i→∞ [g (i) + θg(i)] = 0.

From (2.46) and (2.47), we get

g (i) + θg(i) = (i + 1) ∞ ∑ s=1 ( i + 1 i + 1 + s ) θ -1 + θs i + 1 i ∑ w=0 λ t-1 (w) i-w ∑ h=0 B i-w (h) min( ĥ,L t ) ∑ n=0 λ t (n + s)ξ ĥn .
Since ( i+1 i+1+s ) θ -1 + θs i+1 ≥ 0, after some algebraic operations, we have

g (i) + θg(i) ≤ (i + 1) ∞ ∑ n=1 λ t (n) n ∑ s=1 ( i + 1 i + 1 + s ) θ -1 + θs i + 1 • ξ ĥ,n-s
Using the following inequalities for v ≥ 0 and 0

< θ < 1, 0 ≤ 1 (1 + v) θ -1 + θv ≤ θ(1 + θ) v 2 2 , (2.50)
we have

0 ≤g (i) + θg(i) ≤ θ(θ + 1) 2 (i + 1) i ∑ w=0 λ t-1 (w) i-w ∑ h=0 B i-w (h) N ∑ n=1 λ t (n) n ∑ s=1 s 2 (i + 1) 2 ξ ĥ,n-s + θ(i + 1) i ∑ w=0 λ t-1 (w) i-w ∑ h=0 B i-w (h) ∞ ∑ n=N+1 λ t (n) n ∑ s=1 s i + 1 ξ ĥ,n-s ≤ 1 2(i + 1) N ∑ n=1 n 2 λ t (n) + ∞ ∑ n=N+1 nλ t (n). When L t =o((lni) 4 ) * , according to (2.23) in Lemma 2.4, we can choose N = i 1/3 L t Λ such that P{N t ≥ i 1/3 } ≤ e -λ (eλ) i 1/3 (i 1/3 ) i 1/3 ≤ e -i 1/3 ( L t Λ i 1/3 -1) L t Λ i 1/3 i 1/3 ≤ 1 a i 1/3 4 , (2.51)
where a 4

i 1/3 ΛL t e ΛL t i 1/3 -1 > 1 2 i 1/3 . Then, fix N = i 1/3 , we have 1 2(i + 1) N ∑ n=1 n 2 λ t (n) ≤ i 1/3 2(i + 1) L t Λ ≤ i -1/3 . (2.52)
As a consequence, we have lim i→∞ [g (i) + θg(i)] = 0 at last, which completes the proof of the second part and also Lemma 2.6.

Proof of Lemma 2.7

Proof. Recall R1 , we have R1 ≤ ∑ M x 0 =1 x 0 e -α α x 0

x 0 ! Φ(α). We then rewrite the upper bound Φ(α) of R1 as

Φ(α) = e -α M ∑ m=1 α m (m -1)! ,
whose derivative can be calculated as

Φ (α) = e -α M-1 ∑ m=0 α m m! - α M (M -1)! . ( 2

.53)

We distinguish two cases to look for α * .

Case 1: α ≥ M.

Since it holds that N! ≤ N N-1 for ∀N ∈ N, we can get

Φ (α) < e -α (M -1)! • M ∑ m=1 M M-m α m-1 -α M < e -α (M -1)! • Mα M-1 -α MM < 0,
meaning that Φ(α) monotonously decreases when α≥M.

Case 2: α ≤ M-1 e . Substituting the inequality N! ≥ ( N e ) N into (2.53) yields

Φ (α) ≥ e -α (M -1)! • ( M -1 e ) M-1 + α( M -1 e ) M-1 -α M ≥ e -α (M -1)! • α M-1 + α M -α M > 0,
meaning that Φ(α) monotonously increases as α≤ M-1 e . Combining the analysis in both cases, we have proved that α * maximising Φ(α) falls into the interval

( M-1 e , M), i.e., α * = Θ(M).

Context and Motivation

Recent years have witnessed an unprecedented development and application of the radio frequency identification (RFID) technology. As a promising low-cost technology, RFID is widely utilized in various applications ranging from inventory control [START_REF][END_REF][39], supply chain management [START_REF] Lee | Efficient storage scheme and query processing for supply chain management using RFID[END_REF] to tracking/location [41][42]. A standard RFID system has two types of devices: a set of RFID tags and one or multiple RFID readers (simply called tags and readers). A tag is typically a low-cost microchip labeled with a unique serial number (ID) to identify an object. A reader, on the other hand, is equipped with an antenna and can collect the information of tags within its coverage area.

Tag population estimation and counting is a fundamental functionality for many RFID applications such as warehouse management, inventory control and tag identification. For example, quickly and accurately estimating the number of tagged objects is crucial in establishing inventory reports for large retailers such as Wal-Mart [START_REF]Wal-Mart begins RFID process changes[END_REF]. Due to the paramount practical importance of tag population estimation, a large body of studies [START_REF] Kodialam | Anonymous tracking using RFID tags[END_REF] [45] [START_REF] Qian | Cardinality estimation for large-scale RFID systems[END_REF] [47] [START_REF] Zheng | Zoe: Fast cardinality estimation for large-scale RFID systems[END_REF] have been devoted to the design of efficient estimation algorithms. Most of them, as reviewed in Sec. 3.2, are focused on the static scenario where the tag population is constant during the estimation process. However, many practical RFID applications, such as logistic control, are dynamic in the sense that tags may be activated or terminated as specialized in C1G2 standard [START_REF] Inc | Radio-frequency identity protocols class-1 generation-2 UHF RFID protocol for communications at 860 mhz -960 mhz version 1[END_REF], or the tagged objects may enter and/or leave the reader's covered area frequently, thus resulting in tag population variation. In such dynamic applications, a fundamental research question is how to design efficient algorithms to dynamically trace the tag

Summary of Contributions

In this chapter, we develop a generic framework of stable and accurate tag population estimation schemes for both static and dynamic RFID systems. By generic, we mean that our framework both supports the real-time monitoring and can estimate the number of tags accurately without any prior knowledge on the tag arrival and departure patterns. Our design is based on the extended Kalman filter (EKF) [START_REF] Song | The extended Kalman filter as a local asymptotic observer for nonlinear discretetime systems[END_REF], a powerful tool in optimal estimation and system control. By performing Lyapunov drift analysis, we mathematically prove the efficiency and stability of our framework in terms of the boundedness of estimation error.

The main technical contributions of this chapter are articulated as follows. We formulate the system dynamics of the tag population for both static and dynamic RFID systems where the number of tags remains constant and varies during the estimation process. We design an EKF-based population estimation algorithm for static RFID systems and further enhance it to dynamic RFID systems by leveraging the cumulative sum control chart (CUSUM) to detect the population change. By employing Lyapunov drift analysis, we mathematically characterise the performance of the proposed framework in terms of estimation accuracy and convergence speed by deriving the closed-form conditions on the design parameters under which our scheme can stabilise around the real population size with bounded relative estimation error that tends to zero within exponential convergence rate. To the best of our knowledge, our work is the first theoretical framework that dynamically traces the tag population with closed form conditions on the estimation stability and accuracy.

Related Work

Due to its fundamental importance, tag population estimation has received significant research attention, which we briefly review in this section.

Tag Population Estimation for Static RFID systems

Most of existing works are focused on the static scenario where the tag population is constant during the estimation process. The central question there is to design efficient algorithms quickly and accurately estimating the static tag population. Kodialam et al. design an estimator called PZE which uses the probabilistic properties of empty and collision slots to estimate the tag population size [START_REF] Kodialam | Fast and reliable estimation schemes in RFID systems[END_REF]. The authors then further enhance PZE by taking the average of the probability of idle slots in multiple frames as an estimator in order to eliminate the constant additive bias [START_REF] Kodialam | Anonymous tracking using RFID tags[END_REF]. Han et al. exploit the average number of idle slots before the first non-empty slots to estimate the tag population size [START_REF] Han | Counting RFID tags efficiently and anonymously[END_REF]. Later, Qian et al. develop Lottery-Frame scheme that employs geometrically distributed hash function such that the jth slot is chosen with prob. 1 2 j+1 [START_REF] Qian | Cardinality estimation for large-scale RFID systems[END_REF]. As a result, the first idle slot approaches around the logarithm of the tag population and the frame size can be reduced to the logarithm of the tag population, thus reducing the estimation time. Subsequently, a new estimation scheme called ART is proposed in [START_REF] Shahzad | Every bit counts: fast and scalable RFID estimation[END_REF] based on the average length of consecutive non-empty slots. The design rational of ART is that the average length of consecutive non-empty slots is correlated to the tag population. ART is shown to have smaller variance than prior schemes. More recently, Zheng et al. propose another estimation algorithm, ZOE, where each frame just has a single slot and the random variable indicating whether a slot is idle follows Bernoulli distribution [START_REF] Zheng | Zoe: Fast cardinality estimation for large-scale RFID systems[END_REF]. The average of multiple individual observations is used to estimate the tag population.

We would like to point out that the above research work does not consider the estimation problem for dynamic RFID systems and thus may fail to monitor the system dynamics in real time. Specifically, in typical static tag population estimation schemes, the final estimation result is the average of the outputs of multi-round executions. When applied to dynamic tag population estimation, additional estimation error occurs due to the variation of the tag population size during the estimation process.

Tag Population Estimation for Dynamic RFID systems

Only a few propositions have tackled the dynamic scenario. The works in [START_REF] Sarangan | A framework for fast RFID tag reading in static and mobile environments[END_REF] and [START_REF] Xie | Efficient tag identification in mobile RFID systems[END_REF] consider specific tag mobility patterns that the tags move along the conveyor in a constant speed, while tags may move in and out by different workers from different positions, so these two algorithm cannot be applicable to generic dynamic scenarios. Subsequently, Xiao et al. develop a differential estimation algorithm, ZDE, in dynamic RFID systems to estimate the number of arriving and removed tags [START_REF] Xiao | Differential estimation in dynamic RFID systems[END_REF]. More recently, they further generalize ZDE by taking into account the snapshots of variable frame sizes [START_REF] Xiao | Temporally or spatially dispersed joint rfid estimation using snapshots of variable lengths[END_REF]. Though the algorithms in [START_REF] Xiao | Differential estimation in dynamic RFID systems[END_REF] and [START_REF] Xiao | Temporally or spatially dispersed joint rfid estimation using snapshots of variable lengths[END_REF] can monitor the dynamic RFID systems, they may fail to estimate the tag population size accurately, because they must use the same hash seed in the whole monitoring process, which cannot reduce the estimation variance. Using the same seed is an effective way in tracing tag departure and arrival, however, it may significantly limit the estimation accuracy, even in the static case.

Besides the limitations above, prior works do not provide formal analysis on the stability and the convergence rate. To fill this vide, we develop a generic framework for tag population estimation in dynamic RFID systems. By generic, we mean that our framework can both support real-time monitoring and estimate the number of tags accurately without the requirement for any prior knowledge on the tag arrival and departure patterns. As another distinguished feature, the efficiency and stability of our framework is mathematically established.

Technical Preliminaries

In this section, we briefly introduce the extended Kalman filter and some fundamental concepts and results in stochastic process which are useful in the subsequent analysis. The main notations used in this chapter are listed in Table 3.1.

Extended Kalman Filter

The extended Kalman filter is a powerful tool to estimate system state in nonlinear discrete-time systems.

Formally, a nonlinear discrete-time system can be described as follows:

z k+1 = f (z k , x k ) + w * k (3.1)
y k = h(z k ) + u * k , (3.2)
where z k+1 ∈ R n denotes the state of the system, x k ∈ R d is the controlled inputs and y k ∈ R m stands for the measurement observed from the system. The uncorrelated stochastic variables w * k ∈ R n and u * k ∈ R m denote the process noise and the measurement noise, respectively. The functions f and h are assumed to be the continuously differentiable. For the above system, we introduce an EKF-based state estimator given in Definition 3.1. [START_REF] Song | The extended Kalman filter as a local asymptotic observer for nonlinear discretetime systems[END_REF]). A two-step discrete-time extended Kalman filter consists of state prediction and measurement update, defined as 1) Time update (prediction)

Definition 3.1 (Extended Kalman filter

ẑk+1|k = f ( ẑk|k , x k ) (3.3) P k+1|k = P k|k + Q k , (3.4)
2) Measurement update (correction)

ẑk+1|k+1 = f ( ẑk+1|k , x k ) + K k+1 v k+1 (3.5
)

P k+1|k+1 = P k+1|k (1 -K k+1 C k+1 ) (3.6) K k+1 = P k+1|k C k+1 P k+1|k C k+1 2 + R k+1 , (3.7)
where

v k+1 = y k+1 -h( ẑk+1|k ) (3.8) C k+1 = ∂h(z k+1 ) ∂z k+1 z k+1 = ẑk+1|k . ( 3.9)
Remark. In the above definition of extended Kalman filter, the parameters can be interpreted in our context as follows:

• ẑk+1|k is the prediction of z k+1 at the beginning of frame k + 1 given by the previous state estimate, while ẑk+1|k+1

is the estimate of z k+1 after the adjustment based on the measure at the end of frame k + 1.

• v k+1 , referred to as innovation, is the measurement residual in frame k+1. It represents the estimated error of the measure.

• K k+1 is the Kalman gain. With reference to equation (3.5), it weighs the innovation v k+1 w.r.t. f ( ẑk+1|k , x k ).

• P k+1|k and P k+1|k+1 , in contrast to the linear case, are not equal to the covariance of estimation error of the system state. Here, we will refer to them as pseudo-covariance.

• Q k and R k are two tunable parameters which play the role as that of the covariance of the process and measurement noises in linear stochastic systems to achieve optimal filtering in the maximum likelihood sense. We will show later that Q k and R k also play an important role in improving the stability and convergence of our EKF-based estimators.

Boundedness of Stochastic Process

In order to analyse the stability of an estimation algorithm, we need to check the boundedness of the estimation error defined as follows:

e k|k-1 z k -ẑk|k-1 . (3.10)
Due to probabilistic nature of the estimation algorithm, the estimation process is a stochastic process. Thus, we further introduce the following two mathematical definitions [START_REF] Morozan | Boundedness properties for stochastic systems[END_REF] [58] on the boundedness of stochastic process.

< ζ < 1 such that E[e 2 k|k-1 ] ≤ ψ 1 e 2 1|0 ζ k-1 + ψ 2 . (3.12)
To investigate the boundedness defined in Definition 3.2 and 3.3, we introduce the following lemma [START_REF] Reif | Stochastic stability of the discrete-time extended Kalman filter[END_REF].

Lemma 3.1. Given a stochastic process V k (e k|k-1 ) and constants β, β, τ>0 and 0<α≤1 with the following properties:

βe 2 k|k-1 ≤ V k (e k|k-1 ) ≤ βe 2 k|k-1 , (3.13) E[V k+1 (e k+1|k )|e k|k-1 ] -V k (e k|k-1 ) ≤ -αV k (e k|k-1 ) + τ, (3.14)
then for any k ≥ 1 it holds that

• the stochastic process e k|k-1 is exponentially bounded in the mean square, i.e.,

E[e 2 k|k-1 ] ≤ β β E[e 2 1|0 ](1 -α) k-1 + τ β k-2 ∑ j=1 (1 -α) j ≤ β β E[e 2 1|0 ](1 -α) k-1 + τ βα , (3.15)
• the stochastic process e k|k-1 is bounded w.p.o..

From Lemma 3.1, it can be known that if we can construct V k (e k|k-1 ), a function of e k|k-1 , such that both its drift and

V k (e k|k-1 ) e 2 k|k-1
are bounded, i.e, (3.14) and (3.13) hold, then e k|k-1 is also bounded and the convergence rate depends on constant α mostly. Besides, it can be noted that Lemma 3.1 can only be implemented offline.

To address this limit, we adjust Lemma 3.1 to an online version with time-varying parameters, which can be proven by the same method as in [START_REF] Tarn | Observers for nonlinear stochastic systems[END_REF] and [START_REF] Rhudy | Online stochastic convergence analysis of the Kalman filter[END_REF].

Lemma 3.2. If there exist a stochastic process V k (e k|k-1 ) and real numbers β * , β k , τ k >0 and 0<α * k ≤1 with the following properties:

V 1 (e 1|0 ) ≤ β * e 2 1|0 , (3.16
)

β k e 2 k|k-1 ≤ V k (e k|k-1 ), (3.17
)

E[V k+1 (e k+1|k )|e k|k-1 ] -V k (e k|k-1 ) ≤ -α * k V k (e k|k-1 ) + τ k ; (3.18)
then for any k ≥ 1 it holds that

• the stochastic process e k|k-1 is exponentially bounded in the mean square, i.e.,

E[e 2 k|k-1 ] ≤ β * β k E[e 1|0 2 ] k-1 ∏ i=1 (1 -α * i ) + 1 β k k-2 ∑ i=1 τ k-i-1 i ∏ j=1 (1 -α * k-j ), (3.19)
• the stochastic process e k|k-1 is bounded w.p.o..

Remark

System Model and Problem Formulation

System Model

Consider a RFID system consisting of a reader and a mass of tags operating on one frequency channel. The number of tags is unknown a priori and can be constant or dynamic (time-varying), which we refer to as static and dynamic systems, respectively throughout this chapter. The MAC protocol for the RFID system is the standard framed-slotted ALOHA protocol, where the standard Listen-before-Talk mechanism is employed by the tags to respond the reader's interrogation [START_REF] Finkenzelle | RFID handbook: Radio frequency identification fundamentals and applications[END_REF].

The reader initiates a series of frames indexed by an integer k ∈ Z + . Each individual frame, referred to as a round, consists of a number of slots. The reader starts frame k by broadcasting a begin-round command with frame size L k , persistence probability r k and a random seed Rs k . When a tag receives a begin-round command, it uses a hash function h(•), L k , Rs k , and its ID to generate a random number i in the range [0, L k -1] and reply in slot i of frame k with probability r k . 1Since every tag picks its own response slot individually, there may be zero, one, or more than one tags transmitting in a slot, which are referred to as idle, singleton, and collision slots, respectively. The reader is not assumed to be able to distinguish between a singleton or a collision slot, but it can detect an idle slot. We term both singleton and collision slots as occupied slots throughout this chapter. By collecting all replies in a frame, the reader can generate a bit-string B k illustrated as

B k = {• • • |0|0|1|0|1|1| • • • }, where '0' indicates an idle slot,
and '1' stands for an occupied one.

Subsequently, the reader finalizes the current frame by sending an end round command. Based on the number of idle slots, i.e., the number of '0' in B k , the reader runs the estimation algorithm, detailed in the next section, to trace the tag population.

Tag Population Estimation Problem

Our objective is to design a stable and accurate tag population estimation algorithm for both static and dynamic systems. By stable and accurate we mean that

• the estimation error of our algorithm is bounded in the sense of Definition 3.2 and 3.3 and the relative estimation error tends to zero;

• the estimated population size converges to the real value with exponential rate.

Mathematically, we consider a large-scale RFID system of a reader and a set of tags with the unknown size z k in frame k which can be static or dynamic. Denote by ẑk|k-1 the prior estimate of z k in the beginning of frame k. At the end of frame k, the reader updates the estimate ẑk|k-1 to ẑk|k by running the estimation algorithm. Our designed estimation scheme need to guarantee the following properties:

• lim

z k →∞ ẑk|k-1 -z k z k = 0;
• the converges rate is exponential.

Tag Population Estimation: Static Systems

In this section, we focus on the baseline scenario of static systems where the tag population is constant during the estimation process. We first establish the discrete-time model for the system dynamics and the measurement model using the bit string B k observed during frame k. We then present our EKF-based estimation algorithm.

System Dynamics and Measurement Model

Consider the static RFID systems where the tag population stays constant, the system state evolves as

z k+1 = z k , (3.20)
meaning that the number of tags z k+1 in the system in frame k + 1 equals that in frame k.

In order to estimate z k , we leverage the measurement on the number of idle slots during a frame. To start, we study the stochastic characteristics of the number of idle slots.

Assume that the initial tag population z 0 falls in the interval z 0 ∈ [z 0 , z 0 ], yet the exact value of z 0 is unknown and should be estimated. The range [z 0 , z 0 ] can be a very coarse estimation that can be obtained by any existing population estimation method. Recall the system model that in frame k, the reader probes the tags with the frame size L k . Denote by variable N k the number of idle slots in frame k, that is, the number of '0's in B k , we have the following results on N k according to [START_REF] Kodialam | Fast and reliable estimation schemes in RFID systems[END_REF][START_REF] Kolchin | Random allocation[END_REF].

= L k (1 -r k L k ) z k and σ 2 = L k (L k -1)(1 -2r k L k ) z k + L k (1 -r k L k ) z k -L k 2 (1 -r k L k ) 2z k .
Lemma 3.4. For any * > 0, there exists some M > 0,

such that if z k ≥ M or L k = ẑk|k-1 ≥ M, then it holds that µ -L k e -r k ρ ≤ * , (3.21) σ 2 -L k (e -r k ρ -(1 + r 2 k ρ)e -2r k ρ ) ≤ * , (3.22)
where ρ = z k L k is referred to as the reader load factor.

Lemmas 3.3 and 3.4 imply that in large-scale RFID systems, we can use L k e -r k ρ and L k (e -r k ρ -(1 + r 2 k ρ)e -2r k ρ ) to approximate µ and σ 2 .

At the end of each frame k, the reader gets a measure y k of the idle slot frequency defined as

y k = N k L k . ( 3.23)
Recall Lemma 3.3, it holds that y k is a Normal distributed random variable specified as follows:

E[y k ] = e -r k ρ and Var[y k ] = 1 L k (e -r k ρ -(1 + r 2 k ρ)e -2r k ρ ).
Since there are z k tags reply in frame k with probability r k , the probability that a slot is idle, denoted as p(z k ), can be calculated as

p(z k ) = (1 - r k L k ) z k ≈ e - r k z k L k . (3.24)
Notice that for large z k , p(z k ) can be regarded as a continuously differentiable function of z k .

Using the language in the Kalman filter, we can write y k as follows:

y k = p(z k ) + u k , (3.25)
where, based on the statistic characteristics of y k , u k is a Gaussian random variable with zero mean and variance

Var[u k ] = 1 L k (e -r k ρ -(1 + r 2 k ρ)e -2r k ρ ). (3.26)
We note that u k measures the uncertainty of y k .

To summarise, the discrete-time model for static RFID systems is characterized by (3.20) and (3.25).

Tag Population Estimation Algorithm

Noticing that the system state characterised by (3.20) and (3.25) is a discrete-time nonlinear system, we thus leverage the two-step EKF described in Definition 3.1 to estimate the system state. In (3.7), the Kalman gain K k increases with Q k while decreases with R k . As a result, Q k and R k can be used to tune the EKF such that increasing Q k and/or decreasing R k accelerates the convergence rate but leads to larger estimation error. In our design, we set Q k to a constant q > 0 and introduce a parameter φ k as follows to replace R k to facilitate our demonstration:

R k = φ k P k|k-1 C k 2 . (3.27)
It can be noted from (3.7) and (3.27) that K k is monotonously decreasing in φ k , i.e., a small φ k leads to quick convergence with the price of relatively high estimation error. Hence, choosing the appropriate value for φ k consists of striking a balance between the convergence rate and the estimation error. In our work, we take a dynamic approach by setting φ k to a small value φ but satisfying (3.62) at the first few rounds (J rounds) of estimation to allow the system to act quickly since the estimation in the beginning phase can be very coarse.

After that we set φ k to a relatively high value φ to achieve high estimation accuracy. Now, we present our tag population estimation algorithm in Algorithm 1 where P 0|0 , q can be set to some constants straightforward since the performance mostly depends on φ k and k max is the pre-configured time horizon during which the system needs to be monitored. The major procedures can be summarised as:

1. In the beginning of frame k: prediction (line 3). The reader first predicts the number of tags based on the estimation at the end of frame k-1. The predicted value is defined as ẑk|k-1 . Then the reader sets the persistence probability r k following Lemma 3.8 and z k is set to ẑk|k-1 .

2. Line 4-5. The reader launches the Listen-before-talk protocol as introduced in 3. We will theoretically establish the stability and accuracy of the algorithm in Sec. 3.7.

Algorithm 1 Tag population estimation (static cases): executed by the reader Require: z 0 , P 0|0 , q, J, L, φ, φ, maximum number of rounds k max Ensure: Estimated tag population set

S z = { ẑk|k : k ∈ [0, k max ]} 1: Initialisation: ẑ0|0 ← z 0 , Q 0 ← q, S z = { ẑ0|0 } 2: for k = 1 to k max do 3: ẑk|k-1 ← ẑk-1|k-1 , L k ← L, r k ← 1.59L k / ẑk|k-1 , P k|k-1 ← P k-1|k-1 + Q k-1 4:
Generate a new random seed Rs k and broadcast (L k , r k , Rs k ) 5:

Run Listen-before-Talk protocol 6:

Obtain the number of idle slots N k , and compute y k and v k using (3.23) and (3.8) 7:

Q k ← q 8:
if k ≤ J then S z ← S z ∪ { ẑk|k } 16: end for

Tag Population Estimation: Dynamic Systems

In this section, we further tackle the dynamic case where the tag population may vary during the estimation process. The objective for the dynamic systems is to promptly detect the global tag papulation change and accurately estimate the quantity of this change. To that end, we first establish the system model and then present our estimation algorithm.

System Dynamics and Measurement Model

In dynamic RFID systems, we can formulate the system dynamics as

z k+1 = z k + w k , (3.28)
where the tag population z k+1 in frame k+1 consists of two parts: i) the tag population in frame k and ii) a random variable w k which accounts for the stochastic variation of tag population resulting from the tag arrival/departure during frame k. Notice that w k is referred to as process noise in Kalman filters and the appropriate characterisation of w k is crucial in the design of stable Kalman filters, which will be investigated in detail later. Besides, the measurement model is the same as the static case. Hence, the discrete-time model for dynamic RFID systems can be characterized by (3.28) and (3.25).

Tag Population Estimation Algorithm

In the dynamic case, we leverage the two-step EKF to estimate the system state combined with the CUSUM test to further trace the tag population fluctuation.

Algorithm 2 Tag population estimation (unified framework): executed by the reader Require: z 0 , P 0|0 , q, J, L, φ, φ, maximum number of rounds k max Ensure:

Estimation set S z = { ẑk|k : k ∈ [0, k max ]} 1: Initialisation: ẑ0|0 ← z 0 , Q 0 ← q, S z = { ẑ0|0 } 2: for k = 1 to k max do 3: ẑk|k-1 ← ẑk-1|k-1 , L k ← L, r k ← 1.59L k / ẑk|k-1 , P k|k-1 ← P k-1|k-1 + Q k-1 4:
Generate a new seed Rs k and broadcast (L k , Rs k ) and run Listen-before-Talk protocol 5:

Obtain the number of idle slots N k , and compute y k and v k using (3.23) and (3.8) 6: Our main estimation algorithm is illustrated in Algorithm 2. The difference compared to the static scenario is that tag population variation needs to be detected by the CUSUM test presented in Algorithm 3 in the next subsection and the output of Algorithm 3 acts as a feedback to φ k , meaning φ k is no more a constant after the first J rounds as the static case due to the tag population variation. Specifically, if a change on the tag population is detected in frame k, φ k is set to φ to quickly react to the change, otherwise φ k sticks to φ to stabilize the estimation. The overall structure of the estimation algorithm is illustrated in Fig. 3.1. We note that in the case where z k is constant, Algorithm 2 degenerates to Algorithm 1.

Q k ← q 7: if k ≤ J then 8: φ k ← φ 9:

Detecting Tag Population Change: CUSUM Test

The CUSUM Detection Framework. We leverage the CUSUM test to detect the change of tag population and further adjust φ k . CUSUM test is a sequential analysis technique typically used for change detection [START_REF] Gustafsson | Adaptive filtering and change detection[END_REF]. It is shown to be asymptotically optimal in the sense of the minimum detection time subject to a fixed worst-case expected false alarm rate [START_REF] Brodsky | Nonparametric methods in change point problems[END_REF].

In the context of dynamic tag population detection, the reader monitors the innovation process v k = y kp( ẑk|k-1 ). If the number of the tags population is constant, v k equals to u k which is a Gaussian process with zero mean. In contrast, upon the system state changes, i.e., tag population changes, v k drifts away from the zero mean. In our design, we use Φ k as a normalised input to the CUSUM test by normalising v k with its estimated standard variance, specified as follows:

Φ k = v k (P k|k-1 + Q k-1 )C k 2 + Var[u k ] z k = ẑk|k-1 . (3.29)
The reader further updates the CUSUM statistics g + k and g - k as follows:

g + k = max{0, g + k-1 + Φ k -Υ}, (3.30)
g - k = min{0, g - k-1 + Φ k + Υ}, (3.31)
g + k = g - k = 0, if δ = 1, (3.32)
where g + 0 =0 and g - 0 = 0. And Υ≥0, referred to as reference value, is a filter design parameter indicating the sensitivity of the CUSUM test to the fluctuation of Φ k , Moreover, by δ we define an indicator flag indicating tag population change:

δ =      1 if g + k > θ or g - k < -θ, 0 otherwise, (3.33)
where θ > 0 is a pre-specified CUSUM threshold.

The detailed procedure of the change detection is illustrated in Algorithm 3 where ϕ 1 (δ) is used to assign the value to φ k according to whether the system state changes and is shown in (3.37).

Algorithm 3 CUSUM test: executed by the reader in frame k Require: Υ, θ Ensure: φ k 1: Initialisation:

g + 0 ← 0, g - 0 ← 0 2: Compute Φ k using equation (3.29) 3: g + k ← (3.30), g - k ← (3.31) 4: if g + k > θ or g - k < -θ then 5: δ ← 1, φ k ← ϕ 1 (δ), g + k ← 0, g - k ← 0 6: else 7: δ ← 0, φ k ← ϕ 1 (δ) 8: end if 9: Return φ k
Parameter tuning in CUSUM test. The choice of the threshold θ and the drift parameter Υ has a directly impact on the performance of the CUSUM test in terms of detection delay and false alarm rate. Formally, the average running length (ARL) L(µ * ) is used to denote the duration between two actions [START_REF] Basseville | Detection of abrupt changes: theory and application[END_REF]. For a large θ, L(µ * ) can be approximated as2

L(µ * ) =      Θ(θ), if µ * = 0, Θ(θ 2 ), if µ * = 0, (3.34)
where µ * denotes the mean of the process Φ k .

In our context, ARL corresponds to the mean time between two false alarms in the static case and the mean detection delay of the tag population change in the dynamic case. It is easy to see from (3.34) that a higher value of θ leads to lower false alarm rate at the price of longer detection delay. Therefore, the choices of θ and Υ consists of a tradeoff between the false alarm rate and the detection delay.

Recall that Φ k can be approximated to a white noise process, i.e, Φ k ∼ N [µ * , σ * 2 ] with µ * = 0, σ * = 1 if the system state does not change. Generically, as recommended in [START_REF] Spiring | Introduction to statistical quality control[END_REF], setting θ and Υ as follows achieves good ARL from the engineering perspective.

θ = 4σ * , ( 3.35)
Υ = µ * + 0.5σ * . (3.36)
In the CUSUM framework, we set φ k by ϕ 1 (δ) as follows:

ϕ 1 (δ) =      φ, if δ = 1, φ, if δ = 0. (3.37)
The rationale is that once a change on the tag population is detected in frame k, φ k is set to φ to quickly react to the change, while φ k sticks to φ when no system change is detected.

Performance Analysis

In this section, we establish the stability and the accuracy of our estimation algorithms for both static and dynamic cases.

Static Case

Our analysis is composed of two steps. We first derive the estimation error and then establish the stability and the accuracy of Algorithm 1 in terms of the boundedness of estimation error.

Computing Estimation Error. We first approximate the non-linear discrete system by a linear one. To that end, as the function p(z k ) is continuously differentiable at z k = ẑk|k-1 , using the Taylor expansion, we have

p(z k ) = p( ẑk|k-1 ) + C k (z k -ẑk|k-1 ) + χ(z k , ẑk|k-1 ), (3.38)
where

C k = - r k ρ e r k ρ ẑk|k-1 , (3.39
)

χ(z k , ẑk|k-1 ) = ∞ ∑ j=2 1 e r k ρ j! (r k ρ - r k ρz k ẑk|k-1 ) j . (3.40)
Regarding the convergence of χ(z k , ẑk|k-1 ) in (3.40), assume that

z k = a k ẑk|k-1 , (3.41)
we can obtain the following boundedness of the residual for the case |a k -1| < 1 r k ρ :

|χ(z k , ẑk|k-1 )| = (r k ρ) 2 ( ẑk|k-1 -z k ) 2 e r k ρ ẑ2 k|k-1 ∞ ∑ j=0 (r k ρ) j (j + 2)! 1 - z k ẑk|k-1 j ≤ (r k ρ) 2 ( ẑk|k-1 -z k ) 2 2e (r k ρ) ẑ2 k|k-1 [1 -(r k ρ)(1 -z k ẑk|k-1 ) ] ≤ (r k ρ) 2 ( ẑk|k-1 -z k ) 2 2e (r k ρ) a k ẑ2 k|k-1
, where

a k = 1 -(r k ρ)|1 -a k |. (3.42)
Recall the definition of the estimation error in (3.10) and using (3.20), (3.3) and (3.5), we can derive the estimation error e k+1|k as follows:

e k+1|k =z k+1 -ẑk+1|k = z k -ẑk|k = z k -ẑk|k-1 -K k C k (z k -ẑk|k-1 ) + χ(z k , ẑk|k-1 ) + u k =(1 -K k C k )e k|k-1 + s k + m k , (3.43)
where s k and m k are defined as

s k = -K k u k , (3.44
) 1. there are positive numbers q, q, φ and φ such that the bounds on Q k and φ k are satisfied for every k≥0, as in

m k = -K k χ(z k , ẑk|k-1 ). ( 3
q ≤ Q k ≤ q, (3.46) φ ≤ φ k ≤ φ, (3.47)
2. The initialization must follow the rules

P 0|0 > 0, (3.48) |e 1|0 | ≤ (3.49)
with positive real number > 0.

Remark. By referring to the design objective posed in Section 3.4, Theorem 3.1 prove the following properties of our estimation algorithm:

• the estimation error of our algorithm is bounded in mean square and the relative estimation error tends to zero;

• the estimated population size converges to the real value with exponential rate.

Moreover, the conditions in Theorem 3.1 can be interpreted as follows:

1. The inequalities (3.46) and (3.47) can be satisfied by the configuring the correspondent parameters in Algorithm 1, which guarantees the boundedness of the pseudo-covariance P k|k-1 as shown later.

2. The inequality (3.48) consists of establishing positive P k|k-1 for every k ≥ 1.

3. As a sufficient condition for stability, the upper bound may be too stringent. As shown in the simulation results, stability is still ensured even with a relatively large .

Before the proof of Theorem 3.1, we first state several auxiliary lemmas to streamline the proof and show how to apply these lemmas to prove Theorem 3.1 subsequently.

P k|k-1 = P k-1|k-2 (1 -K k-1 C k-1 ) + Q k-1 = P k-1|k-2 1 - P k-1|k-2 C k-1 2 P k-1|k-2 C k-1 2 + R k-1 + Q k-1 . (3.51)
Following the design of R k in (3.27) and by iteration, we further get

P k|k-1 = P k-1|k-2 1 - 1 1 + φ k-1 + Q k-1 = P 1|0 k-1 ∏ i=1 1 - 1 1 + φ i + k-2 ∑ i=0 Q i k-2 ∏ j=i 1 - 1 1 + φ j+1 + Q k-1 .
Since φ k and Q k are controllable parameters, we can set φ k ≤ φ and Q k ≤ q for every k ≥ 0 in Algorithm 1, where φ, q > 0. Consequently, we have

P k|k-1 ≤ P 1|0 1 - 1 1 + φ k-1 + q k-1 ∑ j=1 1 - 1 1 + φ j + Q k-1 ≤ (P 0|0 + Q 0 ) 1 - 1 1 + φ k-1 + qφ + Q k-1 (3.52) Let p k = ((P 0|0 + Q 0 ) 1 -1 1+φ k-1 + qφ + Q k-1 and p k = Q k-1 , we have p k ≤ P k|k-1 ≤ p k . Lemma 3.6. Let α k 1 1+φ k , it holds that (1 -K k C k ) 2 P k+1|k e 2 k|k-1 ≤ (1 -α k ) e 2 k|k-1 P k|k-1 , ∀k ≥ 1. (3.53)
Proof. From (3.51), we have

P k+1|k = P k|k-1 (1 -K k C k ) + Q k ≥ P k|k-1 (1 -K k C k ) . (3.54)
By substituting it into the left-hand side of (3.53) and using the fact that

R k = φ k P k|k-1 C k 2 for every k ≥ 1, we get (1 -K k C k ) 2 P k+1|k e 2 k|k-1 ≤ (1 -K k C k ) 2 P k|k-1 (1 -K k C k ) e 2 k|k-1 ≤ (1 -K k C k ) e k|k-1 2 P k|k-1 ≤ 1 - 1 1 + φ k e 2 k|k-1 P k|k-1 .
We are thus able to prove (3.53).

Lemma 3.7. Let b k r k ρ(4a k φ k + 1 -a k ) 4a 2 k φ k (1 + φ k ) ẑk|k-1 P k|k-1 , it holds that m k [2(1 -K k C k )e k|k-1 + m k ] P k+1|k ≤ b k | ẑk|k-1 -z k | 3 . (3.55)
Proof. From (3.45), we get the following expansion

m k [2(1 -K k C k )e k|k-1 + m k ] P k+1|k = 1 P k+1|k -P k|k-1 C k P k|k-1 C k 2 + R k χ(z k , ẑk|k-1 ) • 2 1 - P k|k-1 C k 2 P k|k-1 C k 2 + R k e k|k-1 - P k|k-1 C k P k|k-1 C k 2 + R k χ(z k , ẑk|k-1 ) . (3.56)
It then follows from (3.39), (3.41) and (3.54) that

m k [2(1 -K k C k )e k|k-1 + m k ] P k+1|k ≤ 1 P k|k-1 (1 -K k C k ) -P k|k-1 C k P k|k-1 C k 2 + R k (r k ρ) 2 ( ẑk|k-1 -z k ) 2 2e r k ρ a k ẑ2 k|k-1 • 2 1 - P k|k-1 C k 2 P k|k-1 C k 2 + R k | ẑk|k-1 -z k | - P k|k-1 C k P k|k-1 C k 2 + R k (r k ρ) 2 ( ẑk|k-1 -z k ) 2 2e 1.59 a k ẑ2 k|k-1 ≤ r k ρ(4a k φ k + 1 -a k ) 4a 2 k φ k (1 + φ k ) ẑk|k-1 P k|k-1 | ẑk|k-1 -z k | 3 .
We are thus able to prove (3.55).

Lemma 3.8. E s k

2

P k+1|k e k|k-1 ≤ 2.46 ẑk|k-1

φ k (1 + φ k )r k P k|k-1 when r k ρ = 1.59.
Proof. From (3.44), we have

E s 2 k P k+1|k e k|k-1 = K 2 k E[u 2 k ] P k+1|k
. With (3.7), (3.26) and (3.54), we have

E s 2 k P k+1|k e k|k-1 ≤ e 2r k ρ ẑk|k-1 φ k (1 + φ k )P k|k-1 ρr 2 k (e -r k ρ -(1 + r 2 k ρ)e -2r k ρ ). Since item E s 2 k P k+1|k
e k|k-1 influences the estimation accuracy, we set the optimal persistence probability to

minimize this item. Denote Λ(r k ) = e 2r k ρ r 2 k (e -r k ρ -(1 + r 2 k ρ)e -2r k ρ ), we have dΛ dr k = (r k ρ -2)e r k ρ + 2 r 3 k .
Since r k ρ > 0 and d((r k ρ-2)e r k ρ +2)

dr k ρ = (r k ρ -1)e r k ρ which is greater zero if r k ρ > 1 and is smaller than zero if r k ρ < 1, and 1) if r k ρ = 1, dΛ dr k < 0; 2) if r k ρ = 0, dΛ dr k = 0; 3) if r k ρ = 2, dΛ dr k > 0, there exists a unique solution r k ρ ∈ (1, 2) for dΛ dr k = 0 such that Λ(r k ) is minimized. Searching in (1, 2)
, we find the optimal r k ρ = 1.59. Therefore, we can obtain that

E s 2 k P k+1|k e k|k-1 ≤ 2.46 ẑk|k-1 φ k (1 + φ k )r k P k|k-1 ξ k , (3.57)
which completes the proof.

Armed with the above auxiliary lemmas, we next prove Theorem 3.1.

Proof of Theorem 3.1. First, we construct the following Lyapunov function to define the stochastic process:

V k (e k|k-1 ) = e 2 k|k-1 P k|k-1 , which satisfies (3.4) and (3.48) as P k|k-1 > 0.

Next, we use Lemma 3.2 to develop the proof. Because it follows from Lemma 3.5 that the properties (3.16) and (3.17) in Lemma 3.2 are satisfied, the main task left is to prove (3.18).

From (3.43), expanding V k+1 (e k+1|k ) leads to

V k+1 (e k+1|k ) = e 2 k+1|k P k+1|k = [(1 -K k C k )e k|k-1 + s k + m k ] 2 P k+1|k = (1 -K k C k ) 2 P k+1|k e 2
k|k-1

+ m k [2(1 -K k C k )e k|k-1 + m k ] P k+1|k + 2s k [(1 -K k C k )e k|k-1 + m k ] P k+1|k + s 2 k P k+1|k .
Furthermore, by Lemmas 3.6, 3.7 and 3.8 and some algebraic operations, we have

E V k+1 (e k+1|k )|e k|k-1 -V k (e k|k-1 ) ≤ -α k V k (e k|k-1 ) + b k |e k|k-1 | 3 + ξ k . (3.58)
To obtain the same formation with (3.18), we further proceed to bound the second term in b k in (3.58) as follows:

b k |e k|k-1 | 3 ≤ ςα k V k (e k|k-1 ), (3.59)
where 0 < ς < 1 is preset controllable parameter. To prove the above inequality, we need to prove |e k|k-1 | ≤

4ςa 2 k φ k ẑk|k-1 1.59(4a k φ k +1.59|a k -1|) . Since |e k|k-1 | = |a k -1| ẑk|k-1 , it suffices to show |a k -1| ≤ 4ςa 2 k φ k 1.59(4a k φ k + 1.59|a k -1|) , (3.60)
which is equivalent to (1 (3.42). With some algebraic operations, we obtain 1)

-4φ k -4φ k ς)a k 2 + (4φ k -2)a k + 1 ≤ 0 because of
1-2φ k -2 √ φ k (φ k +ς) 1-4φ k (1+ς) < a k ≤ 1, if φ k < 1 4(1+ς) ; and 2) 2φ k -1+2 √ φ k (φ k +ς) 4φ k (1+ς)-1 ≤ a k ≤ 1, if φ k > 1 4(1+ς) ; and 3) 1+ς 1+2ς ≤ a k ≤ 1, if φ k = 1 4(1+ς) . Since it holds that 2φ k -1+2 √ φ k (φ k +ς) 4φ k (1+ς)-1 < 1+ς
1+2ς for every ς and

2φ k -1+2 √ φ k (φ k +ς) 4φ k (1+ς)-1
will decrease monotonically to 1 1+ς for a large φ k , we have in the worst case for φ k ≥ 1 4(1+ς) ,

1 + ς 1 + 2ς ≤ a k ≤ 1. (3.61)
It follows from the analysis that if we set

φ k ≥ 1 4(1 + ς) , (3.62)
(3.60) can be satisfied. Moreover, it holds that

|a k -1| ≤ 0.63ς 1 + 2ς . (3.63)
That is,

|e k|k-1 | ≤ k , (3.64)
where k 0.63ς 1+2ς ẑk|k-1 . By setting φ k in (3.62), for |e k|k-1 | ≤ k , we thus have

E V k+1 (e k+1|k )|e k|k-1 -V k (e k|k-1 ) ≤ -(1 -ς)α k V k (e k|k-1 ) + ξ k . (3.65)
Therefore, we are able to apply Lemma 3.2 to prove Theorem 3.1 by setting = 0.63ς 1+2ς ẑ1|0 ,

β * = 1 Q 0 , α * k = (1 -ς)α k , β k = 1 p k and τ k = ξ k .
Remark. Theorem 3.1 also holds in the sense of Lemma 3.1 (the off-line version of Lemma 3.2) by setting the parameters in (3.15)

as β = 1 Q 0 , α = 1-ς 1+φ ≤ α * k , β = (P 0|0 + Q 0 + q(φ + 1) ≥ p k , and τ = Q 0 ẑmax φ(1+φ) ≥ ξ k ,
where ẑmax is the maximum estimate.

We conclude the analysis on the performance of our estimation algorithm for the static case with a more profound investigation on the evolution of the estimation error |e k|k-1 |. More specifically, we can distinguish three regions:

• Region 1:

2.46M ẑk|k-1 φ k (M-1)r k (1-ς) ≤ |e k|k-1 | ≤ k .
By substituting the condition into the right hand side of (3.65), we obtain:

-(1 -ς)α k V k (e k|k-1 ) + ξ k ≤ -(1-ς)α k M V k (e k|k-1
), where M > 1 is a positive constant and can be set beforehand. It then follows that

E V k+1 (e k+1|k )|e k|k-1 ≤ 1 - (1 -ς)α k M V k (e k|k-1
).

Consequently, we can bound E[e 2 k|k-1 ] as:

E[e 2 k|k-1 ] ≤ p k Q 0 E[e 1|0 2 ] k-1 ∏ i=1 (1 -α * i ) (3.66) with α * k = (1-ς)α k M .
It can then be noted that E[e 2 k|k-1 ] → 0 at an exponential rate as k → ∞.

• Region 2:

2.46 ẑk|k-1 (1-ς) . In this case, we have

φ k r k (1-ς) ≤ |e k|k-1 | < 2.46M ẑk|k-1 φ k (M-1)r k
- (1 -ς)α k M V k (e k|k-1 ) < -(1 -ς)α k V k (e k|k-1 ) + ξ k ≤ 0.
It then follows from Lemma 3.2 that

E[e 2 k|k-1 ] ≤ p k Q 0 E[e 1|0 2 ] k-1 ∏ i=1 (1 -α * i ) + p k k-2 ∑ i=1 ξ k-i-1 i ∏ j=1 (1 -α * k-j ). (3.67) Hence, when k → ∞, E[e 2 k|k-1 ] converges at exponential rate to p k ∑ k-2 i=1 ξ k-i-1 • ∏ i j=1 (1 -α * k-j ) ∼ Θ( ẑk|k-1 ),
which is decoupled with the initial estimation error and it thus holds (1-ς) . In this case, we can show that the right hand side of (3.65) is positive, i.e., -(1ς)α k V k (e k|k-1 ) + ξ k > 0. It also follows from Lemma 3.2 that

E[e k|k-1 ] z k = Θ 1 √ z k → 0 when z k → ∞. • Region 3: 0 ≤ |e k|k-1 | < 2.46 ẑk|k-1 φ k r k
E[e 2 k|k-1 ] ≤ p k Q 0 E[e 1|0 2 ] k-1 ∏ i=1 (1 -α * i ) + p k k-2 ∑ i=1 ξ k-i-1 i ∏ j=1 (1 -α * k-j ). (3.68) Hence, when k → ∞, E[e 2 k|k-1 ] converges exponentially to p k ∑ k-2 i=1 ξ k-i-1 • ∏ i j=1 (1 -α * k-j ) ∼ Θ( ẑk|k-1 ),
which is decoupled with the initial estimation error and it thus holds

E[e k|k-1 ] z k ≤ Θ 1 √ z k → 0 when z k → ∞.
Combining the above three regions, we get the following results on the convergence of the expected estimation error E[e k|k-1 ]: (1) if the estimation error is small (Region 3), it will converge to a value smaller than →0 when z k →∞.

Dynamic Case

Our analysis on the stability of Algorithm 2 for the dynamic case is also composed of two steps. First, we derive the estimation error. Second, we establish the stability and the accuracy of Algorithm 2 in terms of the boundedness of estimation error.

We first derive the dynamics of the estimation error as follows:

e k+1|k = (1 -K k C k )e k|k-1 + s k + m k , (3.69)
which differs from the static case (3.43) in s k . In the dynamic case, we have (3.54) and the independence between w k and e k|k-1 , we can derive

s k = w k -K k u k (3.
> 0 such that E 2s k [(1 -K k C k )e k|k-1 + m k ] P k+1|k e k|k-1 ≤ d k |e k|k-1 |E[w k ]. Proof. When E[w k ] ≥ 0, from E[v k ] = 0, (3.41),
E 2s k [(1 -K k C k )e k|k-1 + m k ] P k+1|k e k|k-1 ≤ 2E[w k ] 1 + φ k φ k P k|k-1 φ k |e k|k-1 | 1 + φ k + 1.59|e k|k-1 | 2 2a k (1 + φ k ) ẑk|k-1 ≤ E[w k ] 2a k φ k + (1 -a k ) a k φ k P k|k-1 |e k|k-1 |.
We thus complete the proof by setting

d k = 2a k φ k + (1 -a k ) a k φ k P k|k-1 . ( 3
φ k (1 + φ k )r k P k|k-1 .
Finally, by setting ξ * k as

ξ * k = 1 + φ k φ k P k|k-1 E[w k 2 ] + 2.46 ẑk|k-1 φ k (1 + φ k )r k P k|k-1 , (3.75)
we complete the proof.

Armed with the above lemmas, we next prove Theorem 3.2 by utilizing the same method with the proof of Theorem 3.1.

Proof of Theorem 3.2. Recall (3.44) and (3.70), we notice that the only difference between the estimation errors of Algorithms 2 and 1 is s k . Therefore, it suffices to study the impact of w k on V k (e k|k-1 ).

It follows from Lemmas 3.6, 3.7, 3.8, 3.9 and 3.10 that

E V k+1 (e k+1|k )|e k|k-1 -V k (e k|k-1 ) ≤ -α k V k (e k|k-1 ) + b k |e k|k-1 | 3 + d k |e k|k-1 |E[w k ] + ξ * k .
Furthermore, bounding the second item in b k as (3.59) and given φ k in (3.62), yields

E V k+1 (e k+1|k )|e k|k-1 -V k (e k|k-1 ) ≤ -(1 -ς)α k V k (e k|k-1 ) + d k |e k|k-1 |E[w k ] + ξ * k for |e k|k-1 | ≤ k .
And we can thus prove Theorem 3.2 by setting = 0.63ς 1+2ς ẑ1|0 ,

β * = 1 Q 0 , α * k = (1 -ς)α k , τ k = ξ * k + d k |e k|k-1 |λ k and β k = 1 p k .
We conclude the analysis on the performance of our estimation algorithm for the dynamic case with a more profound investigation on the evolution of the estimation error |e k|k-1 | and derive the explicit formulas for λ k and σ k . More specifically, we can distinguish three regions:

• Region 1:

9.84M ẑk|k-1 φ k (M-1)r k (1-ς) ≤ |e k|k-1 | ≤ k .
In this case, the objective is to achieve

E V k+1 (e k+1|k )|e k|k-1 -V k (e k|k-1 ) ≤ - 1 M (1 -ς)α k V k (e k|k-1 ) (3.76) so that E[e 2 k|k-1 ] is bounded as E[e 2 k|k-1 ] ≤ p k Q 0 E[e 1|0 2 ] k-1 ∏ i=1 (1 -α * i ). (3.77) That is, it should hold that d k |e k|k-1 |E[w k ] + ξ * k ≤ M-1 M (1 -ς)α k V k (e k|k-1
). To that end, we firstly let the following inequalities hold

     d k |e k|k-1 |E[w k ] ≤ M-1 2M (1 -ς)α k V k (e k|k-1 ), ξ * k ≤ M-1 2M (1 -ς)α k V k (e k|k-1
).

(3.78)

Secondly, substituting (3.73), (3.75) into (3.78) leads to

E[w k ] ≤ a k φ k (1 -ς)|e k|k-1 | (1 + φ k ) (2a k φ k + 1 -a k ) , (3.79
)

E[w k 2 ] ≤ φ k (M -1)r k (1 -ς)|e k|k-1 | 2 -4.92M ẑk|k-1 2M(1 + φ k ) 2 . (3.80)
Thirdly, let

φ k (M -1)r k (1 -ς)|e k|k-1 | 2 2M(1 + φ k ) 2 ≥ 4.92 ẑk|k-1 (1 + φ k ) 2 , (3.81)
and we thus have

|e k|k-1 | ≥ 9.84M ẑk|k-1 φ k (M -1)r k (1 -ς) ˜ , (3.82)
E[w k 2 ] ≤ 2.46 ẑk|k-1 (1 + φ k ) 2 σ k . (3.83)
The rational behind can be interpreted as follows: i) the right term of (3.80) cannot be less than zero and ii) there always exists the measurement uncertainty in the system. Consequently, the impact of tag population change plus the measurement uncertainty should equal in order of magnitude that of only measurement uncertainty, which can be achieved by establishing However, since a k and a k are unknown a priori, we thus need to transform the right hand side of (3.79) to a computable form. From (3.61), we get 1 a k -1 ≤ ς 1+ς such that it holds for the right hand side of (3.79) that

E[w k 2 ] ≤ K 2 k E[u k 2 ]
a k φ k (1-ς)|e k|k-1 | 3(1+φ k )(2a k φ k +1-a k ) ≥ φ k (1-ς) ˜ 3(1+φ k )(2φ k + ς 1+ς ) . Finally, let E[w k ] ≤ φ k (1 -ς) ˜ 3(1 + φ k ) 2φ k + ς 1+ς λ k , (3.84)
we can establish (3.77) and thus get that E[e 2 k|k-1 ] → 0 at an exponential rate when k → ∞.

• Region 2:

9.84 ẑk|k-1 (1-ς) . Given ˜ , λ k and σ k as in Region 1, in this case, we have

φ k r k (1-ς) ≤ |e k|k-1 | < 9.84M ẑk|k-1 φ k (M-1)r k
-(1 -ς)α k V k (e k|k-1 ) + d k |e k|k-1 |E[w k ] + ξ * k ≤ 0. It then follows from Lemma 3.2 that E[e 2 k|k-1 ] ≤ p k Q 0 E[e 1|0 2 ] k-1 ∏ i=1 (1 -α * i ) + p k k-2 ∑ i=1 τ k-i-1 i ∏ j=1 (1 -α * k-j ). (3.85) Hence, when k → ∞, E[e 2 k|k-1 ] converges exponentially to p k ∑ k-2 i=1 τ k-i-1 • ∏ i j=1 (1 -α * k-j ) ∼ Θ( ẑk|k-1
) and it thus holds that

E[e k|k-1 ] z k = Θ( 1 √ z k ) → 0 for z k → ∞. • Region 3: 0 ≤ |e k|k-1 | < 9.84 ẑk|k-1 φ k r k (1-ς) .
The circumstances in this region are very complicated due to E[w k ] and E[w k 2 ], we here thus just consider the worst case that

E[w k ] = λ k and E[w k 2 ] = σ k . Consequently, we have -(1 -ς)α k V k (e k|k-1 ) + d k |e k|k-1 |E[w k ] + ξ * k > 0,
and it then follows from Lemma 3.2 that

E[e 2 k|k-1 ] ≤ p k Q 0 E[e 1|0 2 ] k-1 ∏ i=1 (1 -α * i ) + p k k-2 ∑ i=1 τ k-i-1 i ∏ j=1 (1 -α * k-j ). (3.86) Hence, when k → ∞, E[e 2 k|k-1 ] converges at exponential rate to p k ∑ k-2 i=1 τ k-i-1 ∏ i j=1 (1 -α * k-j ) ∼ Θ( ẑk|k-1
), and thus

E[e k|k-1 ] z k ≤ Θ( 1 √ z k ) → 0 for z k → ∞.
Note that for the case that E[w k ] < λ k and E[w k 2 ] < σ k , the range of Region 3 will shrink and the range of Region 2 will largen.

Integrating the above three regions, we can get the similar results on the convergence of the expected estimation error E[e k|k-1 ] as in the static case.

Discussions

This section discusses the application of the proposed algorithm to the unreliable channel and multi-reader scenario.

Error-prone channel. The unreliable channel may corrupt a would-be idle slot into a busy slot and vice versa. We consider the random error model as [START_REF] Chen | An efficient tag search protocol in large-scale RFID systems with noisy channel[END_REF]. Let t 0 and t 1 be the false negative rate that a would-be empty slot turns into a busy slot and the false positive rate, respectively. Each parameter without error rate is marked with a superscript * to define its counterpart with error rate t 0 and t 1 . Then, we have

p * (Z k ) = t 1 + (1 -t 0 -t 1 )p(Z k ) (3.87) Var * [u k ] = (1 -t 0 -t 1 ) 2 Var[u k ], (3.88)
and thus get the new Kalman gain K * k as

K * k = 1 (1 -t 0 -t 1 ) K k . (3.89)
It is noted that the ideal channel condition is equivalent to the special case where t 0 =t 1 =0. When the channel is totally random, i.e., t 0 =t 1 =0.5, the noisy will overwhelm the measurement and estimation. Nevertheless, if t 0 , t 1 ∈ (0, 0.5), updating the analysis with the new equations, we find that Theorem 3.1 and 3.2 still holds under the same conditions, meaning that the communication error can be compensated successfully.

Multi-reader case. In multi-reader scenarios, we leverage the same approach as [START_REF] Shahzad | Expecting the unexpected: Fast and reliable detection of missing RFID tags in the wild[END_REF]. The main idea is that a back-end server can be used to synchronize all readers such that the RFID system with multiple readers operates as the single-reader case. Specially, the back-end server calculates all the parameters and sends them to all readers such that they broadcast the same parameters to the tags. Subsequently, each reader sends its bitmap to the back-end server. Then the back-end server applies OR operator on all bitmaps, which eliminates the impact of the duplicate readings of tags in the overlapped interrogation region.

Numerical Analysis

In this section, we conduct extensive simulations to evaluate the performance of the proposed tag population estimation algorithms by focusing on the relative estimation error denoted as

REE k = z k -ẑk|k-1 z k
. Specifically, we simulate in sequence both static and dynamic RFID systems where the initial tag population are z 0 = 10 4 with the following parameters: q = 0.1, P 0|0 = 1, J = 3, θ = 4 and Υ = 0.5 with reference to (3.35) and (3.36), L = 1500, φ = 0.25 and φ = 100 such that (3.62) always holds. Since the proposed algorithms do not require collision detection, we set a slot to 0.4ms as in the EPCglobal C1G2 standard [START_REF] Inc | Radio-frequency identity protocols class-1 generation-2 UHF RFID protocol for communications at 860 mhz -960 mhz version 1[END_REF]. We will discuss the effect of φ and φ on the performance in next section.

Algorithm Verification

In the subsection, we show the impact of φ and φ on the system performance. To that end, with REE 0 =0.5, we keep z k =10 4 in static scenario while the tag population varies in order of magnitude from ẑk|k-1 to 0.4 ẑk|k-1 in different patterns in dynamic scenario. Specifically, we set φ=100 while varying φ = 0.25, 0.5, 1 in Fig. 3.2, 3.3, and fix φ=0.25 with varying φ=1, 10, 100 in Fig. 3.4, 3.5. As shown in the figures, a smaller φ leads to rapider convergence rate while the bigger φ, the smaller the deviation. Thus, we choose φ=0.25 and φ=100 in the rest of the simulation. Moreover, we make the following observations. First, as derived in Theorem 3.2, the estimation is stable and accurate facing to a relative small population change, i.e., around the order of magnitude ẑk|k-1 .

Second, the proposed scheme also functions nicely even when the estimation error is as high as 0.4 ẑk|k-1 tags as shown in Fig. 3.3 and 3.5. This is due to the CUSUM-based change detection which detects state changes promptly such that a small value is set for φ k , leading to rapid convergence rate.

Algorithm Performance

In this section, we evaluate the performance of the proposed EKF-based estimator, referred to as EEKF here, in comparison with [START_REF] Kodialam | Anonymous tracking using RFID tags[END_REF] in static scenario and with [START_REF] Xiao | Temporally or spatially dispersed joint rfid estimation using snapshots of variable lengths[END_REF] in dynamic scenario.

Static System (z k = 10 4 )

We evaluate the performance by varying initial relative error as

• REE 0 = z 0 -ẑ0|0 z 0 = 0.
8 means a large initial estimation error.

• REE 0 = 0.5 means a medium initial estimation error.

• REE 0 = 0.2 implies a small initial estimation error and satisfies (3.64) with 0.5 ≤ ς < 1.

The purpose of the first two cases is to investigate the effectiveness of the estimation in relative large initial estimation errors while the third case intends to verify the analytical results ẑ0|0 z 0 >0.79 as shown in (3.64). Note that EZB uses the optimal persistence probability [START_REF] Kodialam | Anonymous tracking using RFID tags[END_REF] that needs to know coarse range of tag size. And estimation time of EZB increases with the width of the range. Fig. 3.6 illustrates the estimation processes with different initial estimation errors. As shown in the figures, the estimation ẑk|k-1 converges towards the actual number of tags within very short time in all cases, despite the initial estimation error. It is worth noticing that EZB suffers the significant deviation though it is faster than EEKF. In this subsection, we evaluate the performance of EEKF for dynamic systems by comparing with the start-of-the-art solution JREP [START_REF] Xiao | Temporally or spatially dispersed joint rfid estimation using snapshots of variable lengths[END_REF] in terms of execution time to achieve the required accuracy. To that end, we refer to the simulation setting in [START_REF] Xiao | Temporally or spatially dispersed joint rfid estimation using snapshots of variable lengths[END_REF]. Specifically, the initial estimation error is 10%. The tag population size changes by following the normal distribution with the mean of 10000 and the variance of 2000 2 and the accuracy requirement is 95%. By taking 9 samplings, we obtain the results as listed in Table 3.2. As shown in Table 3.2, EEKF is more time-efficient than JREP. This is because the persistence probability in JREP is set to optimise the power-of-two frame size, which increases the variance of the number of empty slots and leads to the performance degradation. In contrast, EEKF can minimise this variance while promptly detecting the tag population changes.

Dynamic system

Conclusion

In this chapter, we have addressed the problem of tag estimation in dynamic RFID systems and designed a generic framework of stable and accurate tag population estimation schemes based on Kalman filter. Technically, we leveraged the techniques in extended Kalman filter (EKF) and cumulative sum control chart (CUSUM)

to estimate tag population for both static and dynamic systems. By employing Lyapunov drift analysis, we missing expected tags in the presence of a large number of unexpected tags in the RFID systems in a reliable and time-efficient way.

Missing expected tag

Unexpected tag

Present expected tag

Prior art and limitation

Prior related work can be classified into three categories from the perspective of detecting missing tags: missing tag detection protocols, tag identification protocols, and tag estimation protocols.

There are two types of missing tag detection protocols: probabilistic [START_REF] Tan | How to monitor for missing RFID tags[END_REF][START_REF] Luo | Probabilistic missing-tag detection and energy-time tradeoff in largescale RFID systems[END_REF][START_REF] Luo | Missing-tag detection and energy-time tradeoff in large-scale RFID systems with unreliable channels[END_REF][START_REF] Shahzad | Expecting the unexpected: Fast and reliable detection of missing RFID tags in the wild[END_REF] and deterministic [START_REF] Li | Identifying the missing tags in a large RFID system[END_REF][START_REF] Zhang | Fast identification of the missing tags in a large RFID system[END_REF][START_REF] Liu | Completely pinpointing the missing RFID tags in a time-efficient way[END_REF]. The probabilistic protocols find out a missing tag event with a certain required probability if the number of missing tags exceeds a given threshold, thus they are more time-efficient but return weaker results in comparison with the deterministic protocols that report all IDs of the missing tags. Actually, they can be used together such that a probabilistic protocol is executed in the first phase as an alarm that reports the absence of tags and then a deterministic protocol is executed to report IDs of missing tags. Unfortunately, all missing tag detection protocols except RUN [START_REF] Shahzad | Expecting the unexpected: Fast and reliable detection of missing RFID tags in the wild[END_REF] work on the hypothesis of a perfect environment without unexpected tags and thus fail to effectively detect missing tags in the presence of unexpected tags. Although RUN [START_REF] Shahzad | Expecting the unexpected: Fast and reliable detection of missing RFID tags in the wild[END_REF] is tailored for missing tag detection in the RFID systems with unexpected tags, all unexpected tags may always participate in the interrogation, which leads to the significant degradation of the performance when the unexpected tag population size scales.

Tag identification protocols [START_REF] Myung | Adaptive splitting protocols for RFID tag collision arbitration[END_REF][START_REF] Namboodiri | Energy-aware tag anticollision protocols for rfid systems[END_REF][START_REF] La Porta | Anticollision protocols for single-reader RFID systems: temporal analysis and optimization[END_REF][START_REF] Shahzad | Probabilistic optimal tree hopping for RFID identification[END_REF] can identify all tags in the interrogation region. To detect missing tags, tag identification protocols can be executed to obtain the IDs of the tags present in the population and then the missing tags can be found out by comparing the collected IDs with those recorded in the database.

However, they are usually time-consuming [START_REF] Li | Identifying the missing tags in a large RFID system[END_REF] and fail to work when it is not allowed to read the IDs of tags due to privacy concern.

Tag estimation protocols [START_REF] Qian | Cardinality estimation for large-scale RFID systems[END_REF][START_REF] Shahzad | Every bit counts: fast and scalable RFID estimation[END_REF][START_REF] Zheng | Zoe: Fast cardinality estimation for large-scale RFID systems[END_REF][START_REF] Chen | Understanding RFID counting protocols[END_REF] are used to estimate the number of tags in the interrogation region.

If many expected tags are absent in RFID systems without unexpected tags, a missing tag event may be detected by comparing the estimation and the number of expected tags stored in the database. However, the estimation error may be misinterpreted as missing tags and cause detection error, especially when there are only a few missing tags. Moreover, the estimation protocol cannot handle the case with a large number of unexpected tags.

Proposed solution and main contributions

Motivated by the detrimental effects of unexpected tags on the performance of missing tag detection, we devise a reliable and time-efficient protocol named Bloom filter-based missing tag detection protocol (BMTD).

Specifically, BMTD consists of two phases, each consisting of a number of rounds.

• In each round of the first phase, the reader fist constructs a Bloom filter by mapping all the expected tag IDs into it such that each tag has multiple representative bits. Then the constructed Bloom filter is broadcasted to all tags. If at least one representative bit of a tag is '0's, it finds itself unexpected and will not participate in the rest of BMTD. Thus, the number of active unexpected tags is considerably reduced.

• Subsequently, in each round of the second phase, the reader constructs a Bloom filter by aggregating the feedbacks from the remaining tags and uses it to check whether any expected tag is absent from the population.

The major contributions of this chapter can be articulated as follows. First, we propose a new solution for the important and challenging problem of missing tag detection in the presence of a large number of unexpected tags by employing Bloom filter to filter out the unexpected tags and then detect the missing tags. Second, we perform the theoretical analysis for determining the optimal parameters used in BMTD that minimize the detection time and also meet the required reliability. Third, we perform extensive simulations to evaluate the performance of BMTD. The results show that BMTD significantly outperforms the state-of-the-art solutions.

Related Work

Extensive research efforts have been devoted to detecting missing tags by using probabilistic method [START_REF] Tan | How to monitor for missing RFID tags[END_REF][START_REF] Luo | Probabilistic missing-tag detection and energy-time tradeoff in largescale RFID systems[END_REF][START_REF] Luo | Missing-tag detection and energy-time tradeoff in large-scale RFID systems with unreliable channels[END_REF][START_REF] Shahzad | Expecting the unexpected: Fast and reliable detection of missing RFID tags in the wild[END_REF] and deterministic method [START_REF] Li | Identifying the missing tags in a large RFID system[END_REF][START_REF] Zhang | Fast identification of the missing tags in a large RFID system[END_REF][START_REF] Liu | Completely pinpointing the missing RFID tags in a time-efficient way[END_REF]. Next, we briefly review the existing solutions of missing tag detection problem and introduce Bloom filter.

Probabilistic protocols

The objective of probabilistic protocols is to detect a missing tag event with a predefined probability. Tan et al. initiate the study of probabilistic detection and propose a solution called TRP in [START_REF] Tan | How to monitor for missing RFID tags[END_REF]. TRP can detect a missing tag event by comparing the pre-computed slots with those picked by the tags in the population.

Different from our BMTD, TRP does not take into account the negative impact of unexpected tags. Follow-up works [START_REF] Luo | Probabilistic missing-tag detection and energy-time tradeoff in largescale RFID systems[END_REF][START_REF] Luo | Missing-tag detection and energy-time tradeoff in large-scale RFID systems with unreliable channels[END_REF] employ multiple seeds to increase the probability of the singleton slot. Same to TRP, they are required to know all the tags in the population. The latest probabilistic protocol called RUN is proposed in [START_REF] Shahzad | Expecting the unexpected: Fast and reliable detection of missing RFID tags in the wild[END_REF].

The difference with previous works lies in that RUN considers the influence of unexpected tags and can work in the environment with unexpected tags. However, RUN does not eliminate the interference of unexpected tags fundamentally such that the false positive probability does not decrease with respect to the unexpected tag population size, which shackles the detection efficiency especially in the presence of a large number of unexpected tags. In addition, the first frame length is set to the double of the cardinality of the expected tag set in RUN, which is not established by theoretical analysis and leads to the failure of estimation method in RUN when the number of the unexpected tags is far larger than that of the expected tags.

Deterministic protocols

The objective of deterministic protocols is to exactly identify which tags are absent. Li et al. develop a series of protocols in [START_REF] Li | Identifying the missing tags in a large RFID system[END_REF] which intend to reduce the radio collision and identify a tag not in the ID level but in the bit level. Subsequently, Zhang et al. propose another series of determine protocols in [START_REF] Zhang | Fast identification of the missing tags in a large RFID system[END_REF] of which the main idea is to store the bitmap of tag responses in all rounds and compare them to determine the present and absent tags. But how to configure the protocol parameters is not theoretically analyzed. More recently, Liu et al. [START_REF] Liu | Completely pinpointing the missing RFID tags in a time-efficient way[END_REF] enhance the work by reconciling both 2-collision and 3-collision slots and filtering the empty slots such that the time efficiency can be improved. None of existing deterministic protocols, however, have been designed to work in the chaotic environment with unexpected tags.

Bloom Filter

A Bloom filter is a randomized data structure that is originally from database contexts [START_REF] Bloom | Space/time trade-offs in hash coding with allowable errors[END_REF][START_REF] Rottenstreich | The variable-increment counting bloom filter[END_REF] and is used to records the members of a set but has attracted much research attention in networking applications [START_REF] Hao | Fast dynamic multiple-set membership testing using combinatorial bloom filters[END_REF].

Specifically, given a set A = {a 1 , a 2 , • • • , a n }, Bloom filter operates as follows:

Initialization: Let a bit array BF be the Bloom filter and the length of BF be f , we initialize BF with zero array. Then the filter is incrementally built by inserting items of A by setting certain bits of BF to 1.

Insertion:

To insert an arbitrary item a i ∈ A, we first need to feed a i to k independent hash functions

h 1 , h 2 , • • • , h k to retrieve k values: h v (a i ) mod f for 1 ≤ v ≤ k, which directs to k positions in BF.
Insertion of a i is then achieved by setting the bits in these k positions to 1.

Query:

To determine whether an item b belongs to A, we can check if b has been inserted into the Bloom filter BF. Achieving this requires b to be inserted by the same hash functions and then we can check every bit

BF[h v (b) mod f ] for 1 ≤ v ≤ k.
If all of k bits are set to 1, the Bloom filter asserts b ∈ A; otherwise, b / ∈ A.

System Model and Problem Formulation

System model

Consider a large RFID system consisting of a single RFID reader and a large number of RFID tags. The reader broadcasts the commands and collects the feedbacks from the tags. In the RFID system, the tags can be either battery-powered active ones or lightweight passive ones that are energized by radio waves emitted from the reader. In this chapter, we first take account of the single-reader case and then extend the proposed protocol to the multi-reader case.

The communications between the readers and the tags follow the Listen-before-talk mechanism [START_REF] Han | Counting RFID tags efficiently and anonymously[END_REF]. During the communications, the tag-to-reader transmission rate and the reader-to-tag transmission rate may differ with each other and are subject to the environment. In practice, the former can be either 40kb/s ∼ 640kb/s in the FM0 encoding format or 5kb/s ∼ 320kb/s in the modulated subcarrier encoding format, while the later is normally about 26.7kb/s ∼ 128kb/s [START_REF] Inc | Radio-frequency identity protocols class-1 generation-2 UHF RFID protocol for communications at 860 mhz -960 mhz version 1[END_REF].

Problem formulation

In the considered RFID system, we use E to denote the set of IDs of the expected tags which are expected to be present in a population and target tags to be monitored. In the RFID system, we assume that an unknown number of tags, m, out of these |E| tags are missing. Note that | • | stands for the cardinality of a set. Denote by E r the set of IDs of the remaining |E|m tags that are actually present in the population. Let U be the set of IDs of unexpected tags within the interrogation region of the reader which does not need to be monitored.

The reader may neither knows exactly the IDs of unexpected tags nor does it know the cardinality of U.

Let M be a threshold on the number of missing expected tags. We use P sys to denote the probability that the reader can detect a missing event. The optimum missing tag detection problem is formally defined as follows. Table 4.1 summaries the main notations used in this chapter.

Bloom Filter-based Missing Tag Detection Protocol

Design rational and protocol overview

To improve the time efficiency of detecting missing tags in the presence of a large number of unexpected tags in the population, we limit the interference of unexpected tags in our protocol. To achieve this goal, we employ a powerful technique called Bloom filter which is a space-efficient probabilistic data structure for representing a set and supporting set membership queries [START_REF] Bloom | Space/time trade-offs in hash coding with allowable errors[END_REF] to rule out the unexpected tags in the set U, which efficiently reduces their interference and thus the overall execution time. Following this idea, we propose a Bloom filterbased Missing Tag Detection protocol (BMTD), by which Bloom filters are sequentially constructed by the reader and by the feedbacks from the active tags in the RFID system. The BMTD consists of two phases: 1) the unexpected tag deactivation phase and 2) the missing tag detection phase.

• The first phase is divided into J rounds where the reader constructs J Bloom filters by mapping the recorded IDs in the reader to deactivate the unexpected tags after identifying them.

• The second phase is divided into W rounds. The reader constructs W Bloom filters according to the responses of the remaining active tags and uses the Bloom filters to detect any missing event. Our protocol either detects a missing event or reports no missing event if the reader does not detect a missing event after W rounds.

We elaborate the design of the BMTD in the rest of this section.

Phase 1: unexpected tag deactivation

In Phase 1, we use Bloom filters to reduce the number of active unexpected tags. Specifically, in the j-th round of Phase 1 (j = 1, 2, ..., J), the reader first constructs a Bloom filtering vector by mapping the expected tags in set U into an l j -bit array using k j hash functions with random seed s j . Here, we denote the l j -bit Bloom filter vector as BF 1,j (E). How the values of l j , k j are chosen and how J is calculated are analysed in Sec. 4.5 on parameter optimisation.

Then, the reader broadcasts the l j -bit Bloom filtering vector, k j and s j to all tags. Upon receiving BF 1,j (E), k j , and s j , each tag maps its ID to k j bits pseudo-randomly at positions h 1 (ID), h 2 (ID), • • • , h k j (ID), and checks the corresponding positions in BF 1,j (E). If all of k j bits are 1, then the tag regards itself expected by the reader.

If any of k j bits is 0, the tag regards that it is unexpected and then remains silent in the rest of the time.

Let U j denote the set of the remaining active unexpected tags after the j-th round of Phase 1, and let

U j ∩ BF 1,j (E)
denote the set of unexpected tags that pass the membership test of BF 1,j (E). Since the Bloom filter has no false negatives, the set of remaining active tags can be represented as E r ∪ U j-1 ∩ BF 1,j (E).

After J rounds when Phase 1 is terminated, the number of remaining active unexpected tags, termed as

|U r |, is |U J ∩ BF 1,J (E)|.
The present tag population size can be written as |E r ∪ U r |. Subsequently, the reader enters Phase 2.

Phase 2: missing tag detection

In the second phase, we still employ Bloom filter to detect a missing tag event. Note that the parameters that the reader broadcasts in each round in Phase 2 except random seeds are identical. In the w-th round of Phase 2 (w = 1, 2, ..., W), the reader first broadcasts the parameters containing the Bloom filter size f w , the number of hash functions R w , and a new random seed d w . How their values are chosen and how W is calculated are analysed in Sec. 4.5 on parameter optimisation.

After receiving the configuration parameters, each tag in the set E r ∪ U r selects R w slots at the indexes

h v (ID) (1 ≤ v ≤ R w )
in the frame of f w slots and transmits a short response at each of the R w corresponding slots. As a consequence, a Bloom filter is formed in the air by the responses from the remaining active tags. In each round, there are two types of slots: empty slots and nonempty slots.

According to the responses from the tags, the reader encodes an f w -bit Bloom filter as follows: If the i-th slot is empty, the reader sets i-th bit of the f w -bit vector to be '0', otherwise '1'. Consequently, a virtual Bloom filter is constructed using which the reader then performs membership test. Let BF 2,w (E r ∪ U r ) denote the constructed Bloom filter in w-th round.

To perform membership test, the reader uses tag IDs from the expected tag set E. Specifically, for each ID in E, the reader maps it into R w bits at positions h v (ID)

(1 ≤ v ≤ R w ) in BF 2,w (E r ∪ U r
). If all of them are '1's, then the tag is regarded as present. Otherwise, the tag is considered to be missing. If a missing event is detected in w-round, the reader terminates the protocol without executing the remaining rounds. Otherwise, the reader initiates a new round until the protocol runs W rounds. If the reader does not detect a missing event after W rounds, it reports no missing event, i.e., the number of missing tags m is less than the threshold M.

An illustrative example of BMTD

We present an illustrative example to show the execution of BMTD. Consider an RFID system with 4 tags. We assume that the reader needs to monitor tag 1 and tag 2 and thus knows their IDs, i.e., E={ID1, ID2}, but it is not aware of the presence of tag 3 and tag 4, who are unexpected, i.e., U={ID3, ID4}. In the example, tag 2 is missing from the population.

As shown in (1) of Fig. 4.2(a), the reader first constructs a Bloom filter BF 1,j (E) by mapping IDs in E and broadcasts a message containing BF 1,j (E) and the values of k j and l j . Here we assume J = 1, k j = 2 and l j = 6.

After receiving BF 1,j (E), each tag checks if it is an expected tag. As shown in (2) of Fig. 4.2(a), tag 1 finds itself expected due to the fact that both h 1 (ID1) and h 2 (ID1) are equal to 1. However, tag 4 realizes that it is unexpected for h 1 (ID4) = 0 and deactivates itself. Different from tag 4, actually unexpected tag 3 passes the test and will participate in the rest of BMTD.

As depicted in (1) of Fig

Performance optimisation and parameter tuning

In this section, we investigate how the parameters in the BMTD are configured to minimise the execution time while ensuring the performance requirement.

Tuning parameters in Phase 1

According to the property of Bloom filter, false negatives are impossible. The false positive rate of the Bloom filter BF 1,j (E) in the j-th round in Phase 1, defined as P 1,j , can be calculated as follows [START_REF] Bloom | Space/time trade-offs in hash coding with allowable errors[END_REF]:

P 1,j = 1 -1 - 1 l j |E|k j k j ≈ (1 -e -|E|k j /l j ) k j . ( 4.1)
By rearranging (4.1), we can express the Bloom filter size in the j-th round as

l j = -|E|k j ln(1 -P 1 k j 1,j ) . ( 4.2)
The total time spent in this round can thus be calculated as l j * t r , where t r denotes the per bit transmission time from reader to tags.

We denote C j the cost to detect and deactivate an unexpected tag as follows:

C j = l j t r |U|(1 -P 1,j ) = -t r |E|k j |U|(1 -P 1,j ) ln(1 -P 1 k j 1,j ) . (4.3)
From the expression of C j , it can be noted that C j represents the average time consumed to detect and deactive an unexpected tag in the j-th round. In our design we minimize C j so as to achieve the optimal time-efficiency.

To minimize C j , we first compute the derivative of C j with respect to k j as follows:

dC j dk j = |E|t r P 1 k j 1,j ln P 1,j -k j (1 -P 1 k j 1,j ) ln(1 -P 1 k j 1,j ) |U|(1 -P 1,j )k j (1 -P 1 k j 1,j ) ln 2 (1 -P 1 k j 1,j ) . ( 4.4)
Furthermore, let dC j dk j = 0, we can obtain

P 1 k j 1,j = 1 2 , (4.5)
and the unique minimiser k * j =

ln P 1,j ln 2 as dC j dk j > 0 when k j > ln p 1,j ln 2 , and dC j dk j < 0 when k j < ln p 1,j ln 2 . Therefore, C j reaches the minimum value when P

1 k * j 1,j = 1 2 .
The optimum Bloom filter size, denoted as l * j , can be computed as

l * j = |E|k * j ln 2 . ( 4.6)
The time spent in the j-th round can be computed as |E|t r k * j ln 2 . Therefore, the total execution time of Phase 1, denoted as T 1 , can be derived as

T 1 = J ∑ j=1 |E|t r k * j ln 2 . ( 4.7)
k * j (1 ≤ j ≤ J), as well as J, are set with the parameters in Phase 2 to minimize the global execution time, as analyzed in Sec. 4.5.3 and Sec. 4.5.4.

Let N * be the number of tags still active after Phase 1 (i.e., J rounds), it holds that

N * = |E| -m + |U r |, (4.8)
where U r is the set of unexpected tags still active after Phase 1. Recall (4.5), the expectation of N * can be derived as

E[N * ] = |E| -m + |U| J ∏ j=1 P 1,j = |E| -m + |U| 1 2 ∑ J j=1 k * j . (4.9)

Tuning parameters in Phase 2

Similar to Phase 1, the false positive rate of the w-th round in Phase 2, defined as P 2,w , can be calculated as

P 2,w = 1 -1 - 1 f w N * R w R w ≈ (1 -e -N * R w / f w ) R w . (4.10)
Therefore, the Bloom filter size is

f w = -N * R w ln(1 -P 1 Rw

2,w )

. Moreover, the probability that at least one missing tag can be detected in w-th round, denoted as P d,w , can be computed as

P d,w = 1 -P m 2,w . (4.11)
Following the analysis above, the probability P sys that the reader is able to detect a missing event after at most W rounds in Phase 2, can thus be written as

P sys = 1 - W ∏ w=1 (1 -P d,w ) = 1 -P mW 2,w . (4.12)
It follows from the system requirement that

P sys = 1 -P mW 2,w = α. (4.13)
As a result, we can obtain

f w = -N * R w ln(1 -(1 -α) 1 mWRw ) . ( 4.14)
In the following lemma, we derive the optimum frame size of the Bloom filter f w which is broadcast by the reader in each round of Phase 2.

Lemma 4.1. Let y WR w , the optimum Bloom filter frame size, denoted by f * w , that achieves the detection requirement while minimising the execution time of Phase 2, is as follows:

f * w = -N * R w ln(1 -(1 -α) 1 my * ) (4.15)
where y * = ln (1-α) m ln 1 2 .

Proof. Denote by f the total length of all W Bloom filters in the second phase, we thus have

f = W ∑ w=1 f w = -N * WR w ln(1 -(1 -α) 1 mWRw ) . ( 4.16)
It can be checked that f depends on the product of W and R w which is the total number of hash functions used in Phase 2. To minimize the execution time, let y WR w , we first calculate the derivation of f with respect to y as follows:

d f dy = N * (1 -α) 1 my ln(1 -α) my(1 -(1 -α) 1 my ) ln 2 (1 -(1 -α) 1 my ) - N * ln(1 - (1 -α) 1 my )
.

Imposing

d f dy = 0 yields y = ln(1 -α) m ln 1 2 .
Moreover, when y < ln(1-α) m ln 1 2 , it holds that

d f dy < 0; when y > ln(1-α) m ln 1 2
, it holds that d f dy > 0. Therefore, f achieves the minimum at y * = ln (1-α) m ln 1 2 . The minimum of f w , denoted by f * w can be computed by injecting y = y * into (4.14). The proof is thus completed.

Remark. As the reader does not have prior knowledge on m, the number of missing tags, in the design of BMTD, we require that the detection performance requirement to be hold for any m ≥ M. Hence, f * w and y * are as follows:

f * w = -N * R w ln(1 -(1 -α) 1 My * ) , (4.17)
where y

* = ln(1 -α) M ln 1 2 , (4.18)
where we use m = M in N * and y * , which is the hardest case. Since N * = |E|m + |U r |, it can be checked that the detection probability P sys is monotonically increasing and P 2,w is monotonically decreasing with respect to the number of missing tags m, meaning that m = M makes the detection hardest and any greater m will ease the hardness, it is thus reasonable to use m = M in the rest of the analysis, because if the reader can detect a missing tag event with probability α when m = M, it will fulfill the detection with probability P sys > α when m > M.

In addition, since y * is the total number of hash functions used in Phase 2 and at least one round is executed so as to detect a missing event, y * needs to be a positive integer. Therefore, we set y * = ln(1-α)

M ln 1 2
, which guarantees the required detection performance requirement. Note that R w and W can be set as arbitrary positive integers.

Under the optimum parameter setting derived above, we can calculate the time needed to execute W rounds of Phase 2, denoted by T 2 , as follows:

T 2 = -t t N * y * ln(1 -(1 -α) 1 My * ) , (4.19)
where t t is the time needed by the tags to transmit one bit to the reader. T 2 sets an upper-bound on the execution time of Phase 2.

Tuning k * j and J to minimize worst-case execution time

In this subsection, we study how to set k * j and J to minimize the worst-case execution time, which corresponds to the experience of the execution time where no missing event is detected and hence all the W rounds in the second round need to be executed. We denote the worst-case execution time by T. In the following theorem, we derive the minimiser of E[T]. Theorem 4.1. Denote x ∑ J j=1 k * j , x need to be set to x * as follows to minimise the worst-case execution time of the BMTD:

x * =          0 |U| ≤ U 0 ln -tr |E| ln(1-(1-α 1 My * )) t t y * |U| ln 2 2 -ln 2 |U| > U 0 , (4.20)
where U 0

|E|t r ln(1-(1-α) 1 My * ) -t t y * ln 2 2
. That is, in regard to minimise the worst-case execution time, when the number of unexpected tags does not exceed a threshold U 0 , Phase 1 is not executed, otherwise Phase 1 is executed with the parameters k * j and J set to

∑ J j=1 k * j = x * .
Proof. Recall the two phases of BMTD and (4.7), we can derive the expectation of T as follows:

E[T] = T 1 + T 2 = J ∑ j=1 |E|t r k * j ln 2 + -t t y * E[N * ] ln(1 -(1 -α) 1 My * ) = |E|t r ln 2 J ∑ j=1 k * j + -t t y * |E| -M + |U| 1 2 ∑ J j=1 k j ln(1 -(1 -α) 1 My * ) . ( 4.21)
From (4.21), it can be noted that E[T] is a function of x = ∑ J j=1 k * j . We then calculate the optimum x * that minimizes E[T]. To that end, we compute the derivation of E[T] with respect to x:

dE[T] dx = |E|t r ln 2 + t t y * |U| ln 2 ln(1 -(1 -α) 1 My * ) 1 2 x . (4.22) Since 1 2 x ≤ 1, it thus holds for all x ≥ 0 that dE[T] dx ≥ 0 if |E|t r ln 2 + t t y * |U| ln 2 ln(1-(1-α) 1 My * ) ≥ 0, i.e., |U| ≤ |E|t r ln(1 -(1 -α) 1 My * ) -t t y * ln 2 2 = U 0 . (4.23)
It is worth noticing that E[T] is a monotonic nondecreasing function in this case with respect to x, we thus set x = 0 to minimize the execution time, which means that if the number of unexpected tags is smaller than the threshold U 0 , we should remove the Phase 1 and only execute Phase 2.

In contrast, if |U| > U 0 , dE[T] dx can be negative, zero, or positive. Setting

dE[T]

dx = 0, the optimal value of x to minimise E[T], defined as x * , can be calculated as

x * = ln -t r |E| ln(1-(1-α) 1 My * ) t t y * |U| ln 2 2 -ln 2 .
Remark. Since x * represents the total number of hash functions used in Phase 1, it needs to be a non-negative integer.

Therefore, we set x * either to its ceiling or floor integer depending on which one leads to a smaller E[T]. The parameters k * j and J are set such that

∑ J j=1 k * j = x * .

Tuning k * j and J to minimize expected detection time

The parameters derived in Theorem 4.1 establish that the BMTD is able to detect a missing event with probability equal to or greater than the system requirement α after W rounds of Phase 2. However, in many practical scenarios, the missing event may be detected in the round w < W when the algorithm can be terminated. In this subsection, we derive the parameter configuration (i.e., k * j and J) that minimises the expected detection time. To that end, we first calculate the probability that at least one of the missing tags can be detected for the first time in a given slot and use it to formulate the expectation of the missing event detection time.

Lemma 4.2. The probability that a missing tag can be detected in a given slot of Phase 1, denoted by q, is as follows:

q = 1 -(1 -(1 -α) 1 y * M ) M N * • 1 -(1 -α) 1 y * M . ( 4

.24)

A loose lower-bound for q, denoted as q min , can be established as follows:

q min = 1 -( 1 2 ) M |E|-M+|U| (1 -(1 -α) 1 y * M ). (4.25)
Proof. A missing tag can be detected in a given slot only when at least one missing tag is hashed to this slot and no tag in E r ∪ U r selects the same location. Consider the hardest case for detecting a missing tag event, i.e., m = M, the probability that at least one missing tag maps to the given slot can be given by 1

-(1 -1 f * w ) MR w .
The probability that no tag in E r ∪ U r maps to that slot is equal to (1 -

1 f * w ) N * R w .
Consequently, multiplying the former by the later leads to q, i.e.:

q = 1 -(1 - 1 f * w ) MR w • (1 - 1 f * w ) N * R w ≈ (1 -e -MRw f * w ) • e -N * Rw f * w = 1 -(1 -(1 -α) 1 y * M ) M N * • (1 -(1 -α) 1 y * M ).
We then derive the lower-bound q min . To that end, noticing that q is negatively correlated with N * which falls into the range |E| -M, |E| -M + |U| , we have

q ≥ 1 -(1 -(1 -α) 1 y * M ) M |E|-M+|U| • (1 -(1 -α) 1 y * M ).
On the other hand, noticing that y * = ln(1-α)

M ln 1 2 ≥ ln(1-α) M ln 1 2 , we have q ≥ q min = 1 -( 1 2 ) M |E|-M+|U| (1 -(1 - α) 1 y * M ).
After calculating q, we next derive the expected missing event detection time, denoted by E[T D ].

Theorem 4.2. The expected missing event detection time E[T D ] is given by the following equation:

E[T D ] = |E|t r x ln 2 + t t |E|-M+|U| ∑ N * =|E|-M 1 -(1 -q) f -f q(1 -q) f q |U| N * -|E| + M 1 2 x N * -|E|+M 1 - 1 2 x |U|-N * +|E|-M . ( 4.26)
Proof. Recall (4.16), it holds that there are f =

-N * y * ln(1-(1-α) 1 My * )
slots in Phase 2. We next calculate the number of slots before detecting the first missing tag. It is easy to check that the event that in slot z the reader detects the first missing tag happens if no missing tags is detected in the first z -1 slots while at least one missing tag is detected in slot z. Let Z denote the random variable of z, we have

P{Z = z} = (1 -q) z-1 * q, (4.27)
which is geometrically distributed.

We can then compute the expectation of Z, conditioned by N * , as follows:

E[Z|N * ] = f ∑ z=1 z • P{Z = z} = 1 -(1 -q) f -f q(1 -q) f q . ( 4.28)
Moreover, it follows from the analysis of Phase 1 that the probability that an unexpected tag is still active after Phase 1 is ∏ J j=1 P 1,j . On the other hand, since U r represents the ID set of active unknown tags after Phase 1, recall (4.5) and ∑ J j=1 k * j = x, we can compute the probability of having u active unexpected tags after Phase 1 as follows:

P{|U r | = u} = |U| u J ∏ j=1 P 1,j u 1 - J ∏ j=1 P 1,j |U|-u = |U| u 1 2 x u 1 - 1 2 x |U|-u .
It can be noted that |U r | follows the binomial distribution. Recall the relationship between N * and |U r | in (4.7), it holds that

E[Z] = |E|-M+|U| ∑ N * =|E|-M E[Z|N * ] |U| N * -|E| + M • 1 2 x N * -|E|+M • 1 - 1 2 x |U|-N * +|E|-M (4.29)
Therefore, E[T D ] can be derived as For any x > 2x 0 , we then have

E[T D ] = T 1 + E[Z] • t t = |E|t r x ln 2 + E[Z] • t t . ( 4
g(x) > |E|t r x ln 2 > 2|E|t r x 0 ln 2 = |E|t r x 0 ln 2 + t t q min ≥ g(x 0 )
The lemma is thus proved.

BMTD parameter setting: summary

We conclude this section by streamlining the procedure of the parameter setting in the BMTD:

J j=1 k * j = x * or ∑ J j=1 k * j = x * e .
Following the above two steps, we can obtain all parameters in the BMTD. Note that R w and W can be picked arbitrarily as long as R w W = y * is satisfied, if we set R w = 1 and W = y * , then the Phase 2 of BMTD is reduced to RUN [START_REF] Shahzad | Expecting the unexpected: Fast and reliable detection of missing RFID tags in the wild[END_REF], RUN is thus a special case of our proposed BMTD.

Cardinality estimation

In order to execute the BMTD, the reader needs to estimate the number of unexpected tags |U|. In our work, we use the SRC estimator which is designed in [START_REF] Chen | Understanding RFID counting protocols[END_REF] and is the current state-of-the-art solution. Denote by |E|m + |U| the estimated total number of tags in the system, then the cardinality |U| can be approximated

as |U| = |E| -m + |U| -|E| if m << |E|, |U|.
Because the number of bits that set to one in Bloom filter is concentrated tightly around the mean [START_REF] Mitzenmacher | Probability and computing: Randomized algorithms and probabilistic analysis[END_REF] and [START_REF] Hao | Building high accuracy bloom filters using partitioned hashing[END_REF], once the estimation |U| is obtained, we can calculate the expectation of N * according to (4.9) with m = M and use it as the estimator of N * .

The SRC estimator consists of two phases: rough estimation and accurate estimation. It is proven in [START_REF] Chen | Understanding RFID counting protocols[END_REF] that SRC can obtain a rough estimation n which at least equals to 0.5(|E|m + |U|) after its first phase. In the second phase, SRC can achieve that the relative estimation error is not greater than which is referred to as confidence range with the settings as follows: the frame size L est =

65

(1-0.04 ) 2 and the persistence probability p pe = min{1, 1.6L est / n}.

We then analyse the overhead introduced to estimate the cardinality of U. As proven in [START_REF] Chen | Understanding RFID counting protocols[END_REF], the overhead of SRC estimator is at most O( 12 + log log(|U| + |E|)), which is moderate for large-scale RFID systems with large |U| and |E|.

Fast detection of missing event

In our estimation approach, we require that m |E|, |U|. In case where m is close to |E|, |U|, the estimation may not be accurate. Luckily, in this case, we can quickly detect a missing event in the cardinality estimation phase due to large m.

Specifically, we analyze the SRC estimator's capability of detecting missing event under large m by comparing the pre-computed slots with those selected by the present tags. Recall the proof of Lemma 4.2, we can derive the detection probability in any given slot, defined as q pre , as

q pre = 1 -1 - p pe L est m * 1 - p pe L est (U+E-m) . (4.31)
Since the detections in different slots are independent of each other, the probability of detecting at least one missing tag event by the SRC estimator can be calculated as 1 -(1q pre ) L est which is a increasing function of m. , SRC is able to detect at least a missing tag event with probability one when m is not less than 100, 200, 600, which means that a missing event is detected by SRC and the reader does not need to invoke the BMTD. In the other side, in the cases that m is less than 100, 200, 600, it holds that

| |U|

|U| -1|≤ 0.138, 0.132, 0.128, respectively. With reference to the conclusion drawn from the Fig. 4.4, the BMTD can tolerate these levels of estimation error.

Sensibility to estimation error

The estimation algorithm we use inevitably introduces error on |U|, which may have a negative impact on the performance of the BMTD. In order to investigate this impact, we next illustrate the sensitivity of the detection time to the estimation error. |U| -1| = 0.5, the departure from the detection time without estimation is only 3%. Therefore, it can be concluded that BMTD is very robust to the estimation error.

Enforcing detection reliability

Estimation error also has impact on the reliability of the BMTD as P sys is calculated base on the estimated cardinality.

To enforce the detection reliability, we introduce more rounds to execute additional Bloom filters. The scheme works as follows: After receiving the Bloom filtering vector constructed by the active tags in the set E r ∪ U r in each round of Phase 2, the reader first counts the actual number of '1' bits in the filtering vector, defined as s 1 and uses it to compute the actual false positive probability, denoted by P2,w , as follows:

P2,w = s 1 f * w , (4.32)
because an arbitrary unexpected tag maps to a '1' bit with a probability of s 1 out of f * w . Following (4.13), we have the observed protocol reliability, denoted by Psys , as follows:

Psys = 1 -PMW 2,w . (4.33)
If Psys < α, the reader adds one more round in Phase 2 to further detect the missing tag event until Psys ≥ α.

Discussion on multi-reader case

In large-scale RFID systems deployed in a large area, multiple readers are thus deployed to ensure the full coverage for a larger number of tags in the interrogation region. In such scenarios, we leverage the approach proposed in [START_REF] Kodialam | Anonymous tracking using RFID tags[END_REF] and employed in [START_REF] Shahzad | Expecting the unexpected: Fast and reliable detection of missing RFID tags in the wild[END_REF]. The main idea is that a back-end server is used to synchronize all readers such that the RFID system with multiple readers operates as the single-reader case.

Specially, the back-end server calculates all the parameters involved in BMTD and constructs Bloom filter and sends them to all readers such that they broadcast the same parameters and Bloom filter to the tags.

Furthermore, each reader sends its individual Bloom filtering vector back to the back-end server. When the back-end server receives all Bloom filtering vectors, it applies logical OR operator on all received Bloom filtering vectors, which eliminates the impact of the duplicate readings of tags in the overlapped interrogation region.

Consequently, a virtual Bloom filter is constructed by the back-end server.

Comparison between BMTD and RUN Comparison under different number of missing tags

In this subsection, we evaluate the performance of BMTD under different number of missing tags, which stands for the effectiveness and efficiency of BMTD. To that end, we set |E| = 1000, |U| = 30000, m = 1 : 50 : 901, α = 0.9 and 0.99. Moreover, we set the threshold to M = 1.

Actual reliability: BMTD achieves the required reliability for any missing tag population size when there are a large number of unexpected tags in the RFID systems. Fig. 4.5(a) and 4.5(b) illustrate the actual reliability of BMTD and RUN for α = 0.9 and 0.99, respectively. It can be observed that both BMTD and RUN achieve the reliability more than that required by the system.

Detection time:

The BMTD outperforms the RUN considerably in terms of detection time even in the scenario with the small number of unexpected tag. Fig. 4.8(a) and 4.8(b) show the detection time for α = 0.9 and 0.99, respectively. As shown in the figures, BTMD is able to save time especially when more unexpected tags are present in the population. Moreover, the increase in detection time of BTMD is more slow than that of RUN. This is due to the ability of BTMD that it can detect the missing tag event when estimating the |U| and determine whether to execute the unexpected tag deactivation phase following Lemma 4.3, which is exactly ignored in RUN.

Comparison under different values of threshold

In this subsection, we evaluate the performance of BMTD under different thresholds, which represents the tolerability of BMTD. To that end, we set |E| = 1000, |U| = 30000, m = 100, α = 0.9 and 0.99. Moreover, we choose such M = 50 : 50 : 300 that the threshold can be greater or smaller than or equal to the number of

Conclusions

This chapter has investigated an important problem of detecting missing tags in the presence of a large number of unexpected tags in large-scale RFID systems. Specifically, we aim at detecting a missing tag event in a reliable and time-efficient way. This chapter has presented a two-phase Bloom filter-based missing tag detection protocol (BMTD). In the first phase, we employed Bloom filter to screen out and then deactivate the unexpected

Design Rational

To design missing tag detection algorithms in the multiple-region multiple-group case, we leverage a powerful technique called Bloom filter which is a space-efficient probabilistic data structure for representing a set and supporting set membership queries [START_REF] Bloom | Space/time trade-offs in hash coding with allowable errors[END_REF] to detect a missing event. In our design, we explore the following three natural ideas, each corresponding to a proposed missing tag detection protocol detailed in the next three sections.

Baseline approach. To enable missing tag detection in the multiple-region multiple-group case, we let the reader use the same Bloom filter parameters in each region for each group of tags and construct the Bloom filter based on the responses from the tags to perform missing event detection. This approach, termed as B-detect, is a direct application of Bloom filter to solve our problem.

Adaptive approach. In the baseline approach B-detect, the reader uses the same parameters in each region, which may not be optimum in the case when tags are not evenly distributed across regions. Motivated by this observation, we develop an adaptive approach, named AB-detect, which enables the reader to use different parameters based on the number of tags in the interrogation region the reader queries. Specifically, for each region r, the reader executes one query, to which tags of all the groups in the region respond. The reader constructs a Bloom filter B r for each region containing the response and aggregates B r (1 ≤ r ≤ R) to form a virtual Bloom filter B AB , based on which it detects missing event for each group.

Group-wise approach. We go further by developing a group-wise approach, referred to as GAB-detect. In GAB-detect, the reader executes G group-wise queries for each region r. Only tags of group g (1 ≤ g ≤ G)

in the interrogation region respond to the g-th query. The reader then constructs a Bloom filter B GAB gr for each group g and aggregates B GAB gr (1 ≤ r ≤ R) to form a virtual Bloom filter B GAB g * using the technique in AB-detect, based on which it detects missing event for group g.

By sequentially analysing the above three approaches and mathematically comparing their performance in terms of execution time, we gradually iron out an optimum detection protocol that works in practice.

The Baseline approach

In the B-detect design to enable missing tag detection in the multiple-region case, we let the reader use the same parameters in each region and construct the Bloom filter based on the responses from the tags to perform missing event detection. Specifically, B-detect consists of two phases, detailed as below.

Protocol Description

Phase 1: Query and feedback collection. The reader performs a query in each region r with the same parameter setting ( f , k, s), where f is the length of the Bloom filter vector, k is the number of independent hash functions used to construct the Bloom filter vector, and s is the seed of the hash functions which is identical for all groups and regions. How their values are chosen is analysed in Sec.

i) = B 1 (i) ∨ • • • ∨ B R (i).
The reader then performs membership test. For each tag in E, the reader maps its ID into k bits at positions (h

v (ID) mod f ) (1 ≤ v ≤ k).
If all the corresponding bits in B are 1, then the tag is regarded as present. Otherwise, the tag is considered to be missing. The reader reports a missing event in group g if the number of missing tags is at least M g and no missing event otherwise.

Performance Optimisation and Parameter Tuning

The execution time of B-detect, defined as T B in number of slots, can be written as

T B = R(t 1 + f δ) R f δ, ( 5.1)
where t 1 denotes the time for the reader to broadcast the query parameters and δ denotes the slot duration which we normalise to 1 for notation conciseness. In a large RFID system, it holds that f t 1 , so we ignore t 1 .

In this subsection, we derive the optimum value of f that minimizes T B .

It is well-known that there is no false negative in the Bloom filter membership test and the false positive rate P f p for an arbitrary group g can be calculated as follows [START_REF] Bloom | Space/time trade-offs in hash coding with allowable errors[END_REF]:

P f p = 1 -1 - 1 f (|E|-m)k k ≈ (1 -e -(|E|-m)k/ f ) k , ( 5.2)
where m = ∑ G g=1 m g denotes the total number of missing tags in all groups. By rearranging (5.2), we can express the Bloom filter size as

f = -(|E| -m)k ln(1 -P 1 k f p ) . ( 5.3)
The following theorem derives the optimal values of f and k in the sense of minimising the execution time.

Theorem 5.1. The optimum size of the Bloom filter and the optimum number of hash functions in B-detect, denoted by f * and k * respectively, that minimize the execution time while satisfying the detection reliability requirement for each group g regardless of m g , are as follows:

f * = (|E| -M) • k * -ln(1 -X 1 k * g * ) , ( 5.4
)

k * = ln 1 -α 1 M g * g * ln 1 2 , ( 5.5)
where M=∑ G g=1 M g , X g 1α 1 Mg g , and g * =arg min g X g .

Proof. Recall the definition of a missing event in group g that at least M g tags are missing, the probability that a missing event can be detected in group g by the reader, defined as P dg , can be computed as

P dg = m g ∑ i=M g m g i (1 -P f p ) i P m g -i f p , (5.6)
and P dg has the following property for any m g ≥ M g :

P dg = (1 -P f p ) M g m g ∑ i=M g m g i (1 -P f p ) i-M g P m g -i f p = (1 -P f p ) M g m g -M g ∑ j=0 m g j + M g (1 -P f p ) j P m g -M g -j f p ≥ (1 -P f p ) M g m g -M g ∑ j=0 m g -M g j (1 -P f p ) j P m g -M g -j f p ≥ (1 -P f p ) M g , ( 5.7)
where the first inequality holds due to the inequality below

( m g j+M g ) ( m g -M g j ) = M g -1 ∏ i=0 m g -i M g + j -i ≥ 1, ∀j ∈ [0, m g -M g ],
where the equality holds when m g = M g .

Hence, to ensure the system requirement P dg ≥ α g regardless of m g , we must ensure the following inequality:

(

-P f p ) M g ≥ α g , or P f p ≤ (1 -α 1 Mg g ). 1
Moreover, since P f p is monotonically decreasing and thus (1 -P f p ) M g is monotonically increasing with respect to the number of missing tags m g , meaning that m g = M g makes the detection hardest and any m g larger than M g will ease the hardness, we thus consider the case where m g = M g for 1 ≤ g ≤ G to meet the detection reliability regardless of m g .

Injecting (5.8) into (5.3) with m g = M g leads to

f ≥ -(|E| -M)k ln 1 -1 -α 1 Mg g 1 k
, where M=∑ G g=1 M g . For clarity, let X g 1α 1 Mg g . Because f needs to be set such that the required detection reliability for any group is achieved and k is identical for all groups, we have:

f = (|E| -M)k -ln[1 -(min 1≤g≤G X g ) 1 k ] . ( 5.9)
Without loss of generality, let g * = arg min g X g and let the derivative of the right hand side of (5.9) with respect to k be 0, we can derive that

k * = ln min g X g ln 1 2 = ln 1 -α 1 M g * g * ln 1 2 .
Lemma 5.1. There is no false negative in AB-detect.

Proof. It suffices to prove that if a tag is present, it holds that

B AB (h v (a) mod f R ) = 1, 1 ≤ v ≤ k,
where a denotes the ID of the tag.

Without loss of generality, assume that the tag a is located in region r. Consider any v ≤ k, let

h v (a) = x + y f r , x, y ∈ N, x < f r . Let c = f R f r .
By definition of B r , PB r and B AB , we have

B AB (x + y f r ) = PB r (x + y f r ) = B r (x) = 1, (5.10)
for ∀y ∈ N, y < c. On the other hand, we have

h v (a) mod f R = x + y f r mod (c f r ) = x + (y mod c) f r .
It then follows from (5.10) that

B AB (h v (a) mod f R ) = 1.
The proof is thus completed.

Performance Optimisation and Parameter Tuning

In this section, we investigate how to tune the parameters in AB-detect to minimise the execution time while ensuring the reliability requirement of each group. We first formulate the false positive rate for each group g, defined as P f p,g . Recall the construction of B AB in AB-detect, the probability that any bit in B AB is zero is

∏ g r=1 1 -1 f r ∑ G g=1 k g n gr
. The false positive rate for group g can then be derived as

P f p,g = 1 - R ∏ r=1 1 - 1 f r ∑ G g=1 k g n gr k g ≈ 1 -e -∑ R r=1 ∑ G g=1 kg ngr fr k g . ( 5

.11)

The following theorem derives the optimal values of f r and k g that minimize the execution time while ensuring the group-wise reliability requirement.

Theorem 5.2. The optimum Bloom filter vector size for the region r and the number of hash functions for the group g, denoted as f * r and k * g , that minimize the execution time while satisfying the detection reliability requirement for each group g regardless of m g , are as follows:

f * r = ∑ G g=1 k * g n gr • ∑ R r=1 ∑ g g=1 k * g n gr min g Y * g , ( 5
.12)

k * g = ln(1 -α 1 Mg g ) ln 1 2 , ( 5.13)
where Y *

g -ln[1 -(1 -α 1 Mg g ) 1 k * g ].
The minimum execution time under the above setting, defined as T * AB , is:

T * AB = 1 min 1≤g≤G Y * g   R ∑ r=1 G ∑ g=1 k * g n gr   2 .
(5.14)

Proof. By the same analysis as the proof of Theorem 5.1, we need to ensure the following inequality:

P dg ≥ (1 -P f p,g ) M g or P f p,g ≤ (1 -α 1 Mg g ). ( 5
∑ r=1 G ∑ g=1 k g n gr f r ≤ -ln[1 -(1 -α 1 Mg g ) 1 kg ], 1 ≤ g ≤ G.
For clarity, let Without loss of generality, let g m = arg min g Y g . It can be checked that

Y g -ln[1 -(1 -α 1 Mg g ) 1
k g ≥ k g m ln (1 -α 1 Mg g ) ln(1 -α 1 Mg g m ) , 1 ≤ g ≤ G.
(5.16)

Next we derive the execution time of AB-detect, defined as T AB . We can write T AB as

T AB = R • G • t 1 + R ∑ r=1 f r R ∑ r=1 f r ,
where t 1 denotes the time for the reader to broadcast protocol parameters including the group ID for each group. In a large RFID system, it holds that f r t 1 . As RGt 1 is constant, finding the optimum k g and f r is equivalent to solving the following optimisation problem:

Minimize: T AB = R ∑ r=1 f r (5.17) Subject to: R ∑ r=1 G ∑ g=1 k g n gr f r ≤ Y g m . ( 5.18)
The corresponding Lagrange function can be defined as

L( f r , λ) = R ∑ r=1 f r + λ R ∑ r=1 G ∑ g=1 k g n gr f r -Y g m .
Solving ∇ f r ,λ = 0 yields the following optimum for f r :

f * r = ∑ G g=1 k g n gr • ∑ R r=1 ∑ G g=1 k g n gr Y g m .
T AB thus achieves its minimum with respect to f r as below:

T * AB = ∑ R r=1 ∑ G g=1 k g n gr • ∑ R r=1 ∑ G g=1 k g n gr Y g m = 1 Y g m   R ∑ r=1 G ∑ g=1 k g n gr   2 .
It can be checked that T AB is monotonously increasing in k g . Recall (5.16), it holds that T AB achieves its minimum as below when the equality in (5.16) holds:

T * AB = min k gm k g m   ∑ R r=1 ∑ G g=1 ln (1-α 1 Mg g ) ln(1-α 1 Mg gm ) n gr   2 Y g m . ( 5.19)
In the above equation,

  ∑ R r=1 ∑ G g=1 ln (1-α 1 Mg g ) ln(1-α

Remark. As k *

g needs to be an integer and f r a power-multiple of two, they need to be rounded to the smallest integer and power-multiple of two not smaller than themselves.

Performance Comparison: B-detect vs. AB-detect

Theorem 5.3. Given the optimum parameters in both B-detect and AB-detect, the following relationship between the minimum execution time of B-detect T * B and that of AB-detect T * AB holds:

1 R ≤ T * AB T * B ≤ 2.
Proof. Recall (5.4), (5.5), (5.13), (5.14) and Y * g in Theorem 5.2, with some algebraic operations, it can be known thatln(1 -X 1 k * g * ) in (5.4) is equal to min g Y * g and k * ≥ k * g for ∀g. We then have

T * AB ≤ k * min g Y * g   R ∑ r=1 G ∑ g=1 n gr   2 . Let T * AB k * min g Y * g ∑ R r=1 ∑ G
g=1 n gr 2 and further recall (5.1), we have

T * AB T * B = ∑ R r=1 ∑ G g=1 n gr 2 R * ∑ R r=1 ∑ G g=1 n gr . Expanding ∑ R r=1 ∑ G g=1 n gr 2 leads to   R ∑ r=1 G ∑ g=1 n gr   2 = R ∑ r=1 G ∑ g=1 n gr + R-1 ∑ i=1 R ∑ r=i+1 2 G ∑ g=1 n gi • G ∑ g=1 n gr ≤ R ∑ r=1 G ∑ g=1 n gr + R-1 ∑ i=1 R ∑ r=i+1 G ∑ g=1 n gi + G ∑ g=1 n gr ≤ R R ∑ r=1 G ∑ g=1 n gr .
To guarantee that f r is power-multiple of two, we need to at most double it. It thus holds that Theorem 5.3 leads to the following engineering implications.

• In the worst case, AB-detect doubles the execution time compared to B-detect;

• In a large asymmetric system where the number of regions R is large, AB-detect can achieve significant performance gain.

The Group-wise Approach

In AB-detect, the reader constructs one Bloom filter that contains the response bits of tags of all groups in the interrogation region. Mixing responses from tags of different group may cause "interference" among groups and thus may increase the detection time for certain groups. Motivated by this observation, we develop a group-wise approach, termed as GAB-detect, in which the reader queries one group each time and constructs group-wise Bloom filters to eliminate the inter-group interference.

Protocol Description

Phase 1: Query and feedback collection. The reader performs G queries in each region r. In the g-th query

(1 ≤ g ≤ G), the reader broadcasts a tetrad (g, k g , f gr , s) where g is the group ID of group g, k g is the number of hash functions used by group g tags, f gr is the Bloom filter size used in region r for group g, s is the hash seed which is identical for all regions and groups. Again, we require f gr to be a power-multiple of two. Without loss of generality, we assume that f g1 ≤ f g2

(1 ≤ r ≤ R -1) to a virtual Bloom filter B GAB g *
for each group g by using the expansion and combination technique in AB-detect. The reader then performs membership test for each group g by using B GAB g * .

Performance Optimisation and Parameter Tuning

In this section, we investigate how to tune protocol parameters in GAB-detect to minimise the execution time while ensuring the reliability requirement of each group. We first derive the false positive rate of GAB-detect for any group g, defined as P f p,g . Recall the construction of B GAB g * , the probability that any bit in B GAB g * is zero is

∏ R r=1 1 -1 f gr k g n gr
. Hence, the false positive rate for group g can be derived as

P f p,g = 1 - R ∏ r=1 1 - 1 f gr k g n gr k g ≈ 1 -e -∑ R r=1 kg ngr fgr k g . ( 5

.20)

The following theorem derives the optimal values of f gr and k g that minimize the execution time while ensuring the group-wise reliability requirement.

Theorem 5.4. The optimum Bloom filter vector size and number of hash functions for group g in region r, denoted as f * gr and k * g , that minimize the execution time while satisfying the detection reliability requirement for each group g regardless of m g , are:

f * gr = √ n gr • ∑ R r=1 √ n gr Z * g , (5.21)
k * g = ln(1 -α 1 Mg g ) ln 1 2 , ( 5.22)
The minimum execution time under the above setting, defined as T * GAB , is:

T * GAB = G ∑ g=1 ∑ R r=1 √ n gr 2 Z * g , (5.23)
where Z

* g ln[1-(1-α 1 Mg g ) 1 k * g ] -k * g .
Proof. By the same analysis as the proof of Theorem 5.1, we need to ensure the following inequality:

P f p,g ≤ (1 -α 1 Mg g ). (5.24)
Injecting (5.20) into (5.24) leads to

R ∑ r=1 k g n gr f gr ≤ -ln[1 -(1 -α 1 Mg g ) 1 kg ] k g . For clarity, let Z g -ln[1-(1-α 1 Mg g ) 1 kg ] k g .
Furthermore, the execution time of GAB-detect, defined as T GAB , can be derived as follows

T GAB = RCt 1 + G ∑ g=1 R ∑ r=1 f gr G ∑ g=1 R ∑ r=1 f gr .
Finding the optimum f gr and k g is equivalent to solving the following optimisation problem: . Injecting k * g into T GAB yields the optimum of T GAB and completes the proof.

Minimize: T GAB = G ∑ g=1 R ∑ r=1 f gr (5.25) Subject to: R ∑ r=1 n gr f gr ≤ Z g , 1 ≤ g ≤ G. ( 5

Performance Comparison: AB-detect vs. GAB-detect

In this section, we compare the execution time of AB-detect and GAB-detect. It follows from x gi x wr + x gr x wi ≥ 2 √ x gi x gr • x wi x wr that φ 2 ir ≤ β 2 ir . We then have

(T * AB -T * GAB ) min g Y * g = 2 R-1 ∑ i=1 R ∑ r=i+1 (β ir -φ ir ) ≥ 0.
The proof is thus completed.

Discussion

In this section, we discuss some implementation issues of our proposed missing tag detection algorithms.

Estimating Tag Population

In our algorithms, the reader needs to estimate the number of tags in n gr in each region and for each group.

This may lead to extra overhead prior to missing tag detection. However, this overhead can be limited as the estimation can be achieved in O(log n gr ) time using state-of-the-art estimation approaches. Specifically, we can apply two types of methods to estimate n gr : single-group estimator and multi-group estimator. In the singlegroup estimator, when staying at region r the reader queries with the group ID g and only the tags from g respond. Then it operates like a single-group system. n gr can be estimated by the methods in [START_REF] Chen | Understanding RFID counting protocols[END_REF]. On the other hand, multi-group estimator estimates multiple group sizes simultaneously by employing the maximum likelihood estimation method as in [START_REF] Luo | An efficient protocol for rfid multigroup threshold-based classification[END_REF], which is time-efficient.

Despite the extra overhead due to estimation of n gr , this estimation phase enables the pre-detection of missing tags if the number of missing tags is important (e.g., due to unexpected loss or accidents). More specifically, the reader can achieve pre-detection by comparing the bitmaps constructed by the tag feedbacks and computed a priori by the reader. If a bit that is 1 in the pre-calculated bitmap by reader but turns out to be 0 in the bitmap of the feedbacks, the reader can identify the absence of tags mapped into this slot. If the number of missing tag for a given group exceeds the threshold, a missing event is reported for the group.

Consequently, the reader may not need to execute the fine-grained detection algorithms as developed in the last three sections since missing tag events have already been detected in the estimation phase, thus reducing the time cost.

We may wonder whether existing tag estimation algorithms can be used to detect the missing tag event. case 1 and case 2, and α g = 0.44 + 0.05(g -1), i.e., 0.449≤α g ≤0.999 in case 3 and case 4. The total number of tags in each region is 12000 and the group size is 12000/G in symmetric scenario. In the asymmetric scenario, on the other hand, the total number of tags is randomly chosen from [1000, 5000] in each of the first R/2 regions and [10000,20000] in the remaining regions, and the group size in the same region is identical. The simulation results are obtained by taking the average of 100 independent trials.

Performance Evaluation

Performance under symmetric scenario Fig. 5.2 depicts the execution time of three protocols under different threshold for the four cases in the symmetric scenario. As shown in the results, globally GAB-detect achieves the best time efficiency and AB-detect outperforms B-detect, especially when the detection reliability for each group varies more significantly, i.e., G = 12. This can be explained as follows: The frame size in AB-detec and B-detec are set based on min g Y * g in Theorem 5.1 and 5.2, which overkills the groups with larger Y * g . In contrast, GAB-detect addresses this limit by eliminating the inter-group interference. We can also observe that in some cases, GAB-detect has longer execution time than AB-detect. This is due to the design requirement that the Bloom filter size needs to be the power-multiple of two. However, globally speaking, GAB-detect still outperforms B-detect. Furthermore, we investigate the actual reliability of the proposed schemes. The results demonstrate that all proposed schemes can detect the missing event with probability one. It can be seen from the four subfigures in Fig. 5.3 that GAB-detect outperforms AB-detect and saves execution time up to 70% in comparison to B-detect. This can be interpreted as follows: In the asymmetric scenarios, the performance gap between AB-detect and B-detect is more significant compared to the symmetric scenario because the frame size in B-detect is identical across the regions regardless of the tag size in an individual region while AB-detect distinguishes the regions with different tag sizes when setting the frame size. Furthermore, we investigate the actual reliability of the proposed schemes and the results show that all proposed schemes can detect the missing event with probability one.

To further evaluate the performance and evaluate the analytical results, we conduct a set of numerical analysis in a even more asymmetric scenario where the tag size is randomly chosen from [START_REF] Song | The extended Kalman filter as a local asymptotic observer for nonlinear discretetime systems[END_REF]100] in each of the first R -1 regions and from [5000, 10000] in the remaining region. As shown in the four subfigures in Fig. 5.4, the performance gain of GAB-detect and AB-detect over B-detect is more remarkable. Specifically, the detection time of B-detect is up to 12.6 times as much as that of GAB-detect and AB-detect.

Impact of nonidentical M g

To comprehensively evaluate the performance, we conduct more numerical analysis in both symmetric and asymmetric scenarios which are same with the previous settings except that R is fixed to 6 and M g = g for In what follows, we discuss a number of open questions we judge pertinent to our work and outline several important potential directions for future research.

Open Questions and Future Work

In this section, we develop the discussion on open issues and questions and future work at three levels: the first level focuses on the new RFID application scenarios, the second one takes the inspiration of the problems we studied to the design of more practical Bloom filter, and the third one focuses on the generalization of our work in a broader context.

Algorithm Design for RFID System With Blocker tags

With the wide adoption of RFID technology, privacy problem attracts more concerns, because RFID tags blindly respond to the interrogation of any RFID reader, even an unauthorized one. An effective solution to this privacy issue is to deploy blocker tags. By programming a set of known RFID tag IDs in a blocker tag, it behave exactly the same as the blocked genuine tags such that the privacy of genuine tags can be protected as the block tag always responds together with its blocking genuine tags which leads to the collision and makes the privacyretrieve attack miscarried.

Although effective on the protection of privacy, the blocker tag also blocks the legal reader from obtaining the information of genuine tags, from which multiple problems arise. Specifically, Consider a large-scale dynamic system where tags may join or leave, a number of blocker tags are deployed to protect the partial tags (e.g., the tags attacked to the most expensive goods) which are referred to as key tags and the other tags are called ordinary tags. In such a dynamic system, the introduction of blocker tags challenges the algorithm design, because blocker tags will mislead the reader in believing the presence of all key tags. Consequently, research effort should be devoted to address the following problems:

• Estimating key tags: a natural and primary concern is how many key tags are present at the current moment.

• Detecting missing key tags: whether there are key tags absent from the population.

• Reconfiguring blocker tags: If the key tags blocked by some blocker tags are not in the systems any more, then how to quickly reconfigure blocker tags to protect the remaining key tags.

Towards Practical Bloom Filter

Our works presented in the thesis can also motivate and inspire a number of design perspectives for Bloom filter.

Fault-tolerable Bloom filter. In Chapter 3, we consider the influence of error-prone channel on estimation efficiency where a would-be empty slot turns into a busy slot and a would-be busy slot turns into an empty slot due to the unreliable channel conditions. While such cases also happen in the applications based on Bloom filter. For example, due to the error in the hardware or the malicious attacks, bit 1 in Bloom filter may become into 0 and vice versa. In such context, it is of great importance to design fault-tolerable Bloom filter, however, the existing work ignores the fault tolerance of Bloom filter.

To address this problem, one direct way is to take into consideration the error rate when configuring the Bloom filter parameters, which, however, is at the price of enlarging the needed store space. Therefore, how to effectively and efficiently design fault-tolerable Bloom filter remains an open research problem.

Heterogeneous Multi-Set Bloom filter. In Chapter 5, we consider a realistic RFID system of multiple tag groups with asymmetric requirement on detection reliability. Such heterogeneity also exists in Bloom filterbased networking applications where data streams of some categories are more popular than the others and thus will be queried more frequently and should be identified with lower false positive rate. The prior work on multi-set Bloom filter, however, does not differentiate this heterogeneity among different categories of data flows. Consequently, how to construct multi-set Bloom filter with the consideration of heterogeneity should be devoted more research effort.

Extension to Big Network Data

In this thesis, we mainly focused on the tag counting and monitoring in RFID systems. A natural extension for our work is to investigate the cardinality estimation (i.e., flow counting and monitoring) in the networks with big data. Currently, the cardinality and categories of data flow are increasing explosively. In this context, it is of great importance to design efficient cardinality estimation algorithm to quickly count the number of flows and data packets in a flow, for practical applications, such as the statistics of popularity of searched keywords in Google and optimisation of traffic routing at ISP side. In a broader context, an important research avenue is to develop efficient algorithms for measuring, monitoring, and analyzing big data transferred over the emerging big data networks such as Internet, social networks or other types of networks. étendu et la technique CUSUM pour estimer le nombre d'étiquettes pour les deux systèmes statiques et dynamiques. En employant l'analyse de Lyapunov, nous caractérisons mathématiquement la performance de notre cadre en termes de la précision de l'estimation et de la vitesse de convergence en dérivant les conditions mathématiques sur les paramètres dans lesquelles notre système peut se stabiliser autour du vrai nombre avec l'erreur relative qui tend vers zéro avec un taux de convergence exponentiel.

Nous procédons en suite à résoudre le problème de détection des étiquettes manquantes, l'une des applica-

Figure 1 . 1 :

 11 Figure 1.1: Thesis organization

  Symposium on Quality of Service (IWQoS), Portland, OR, USA, June 2015. The extended version is currently under second round (minor) revision of IEEE Transactions on Mobile Computing (TMC). Our work on tag population estimation in dynamic RFID systems is currently under second round (minor) revision of IEEE Transactions on Communications (TCOM). Our work on missing tag event detection in the presence of unexpected tags is currently under submission of IEEE Transactions on Communications (TCOM). Our work on missing tag event detection in multiple-group multiple-region RFID systems is under second round (minor) revision of IEEE Transactions on Mobile Computing (TMC).

  )

Figure 2 . 1 :

 21 Figure 2.1: The evolution of backlog population in FSA-SPR.

Figure 2 . 2 :

 22 Figure 2.2: The evolution of backlog population in FSA-MPR.

Figure 2 . 3 :

 23 Figure 2.3: SPR: # of backlogs.

Figure 2 . 4 :

 24 Figure 2.4: MPR: # of backlogs.

Figure 2 . 5 :

 25 Figure 2.5: SPR vs. MPR.

Chapter 3 .

 3 From Static to Dynamic Tag Population Estimation: An Extended Kalman Filter Perspective population quickly and accurately.

Lemma 3 . 3 .

 33 If each tag replies in a random slot among the L k slots with probability r k , then it holds that N k ∼ N [µ, σ 2 ] for large L k and z k , where µ

4 .1 in order to receive the feedbacks from tags. 3 .

 43 At the end of frame k: correction (line 6-14). The reader computes N k based on B k and further calculates y k and v k from N k . It then updates the prediction with the corrected estimate ẑk|k following (3.5).

13 :Figure 3 . 1 :

 1331 Figure 3.1: Estimation algorithm diagram: Dashed box indicates the EKF.

. 45 )Theorem 3 . 1 .

 4531 Boundedness of Estimation Error.Having derived the dynamics of the estimation error, we now state the main result on the stochastic stability and accuracy of Algorithm 1. Consider the discrete-time stochastic system given by (3.20) and (3.25) and Algorithm 1, the estimation error e k|k-1 defined by (3.10) is exponentially bounded in mean square and bounded w.p.o., if the following conditions hold:

  and (3.81) with reference to (3.74) and (3.75).

Figure 3 . 2 :

 32 Figure 3.2: Static:φ = 100.

Figure 3 . 3 :Figure 3 . 4 :Figure 3 . 5 :

 333435 Figure 3.3: Dynamic:φ = 100.

8 Figure 3 . 6 :

 836 Figure 3.6: Algorithm performance under different initial estimation errors.

Figure 4 . 1 :

 41 Figure 4.1: Missing tag detection with the presence of unexpected tags.

Definition 4 . 1 (

 41 Optimum missing tag detection problem). Given |U| unexpected tags where both |U| and the IDs of tags in U are unknown, the optimum missing tag detection problem is to devise a protocol of minimum execution time capable of detecting a missing event with probability P sys ≥ α if m ≥ M, where α is the system requirement on the detection reliability.

Table 4 . 1 :

 41 Main Notations Symbols Descriptions E set of target tags that need to be monitored E r tags that are actually present in the population U set of unexpected tags α required detection reliability m number of missing expected tags M threshold to detect missing tags P sys prob. of detecting a missing event in BMTD J number of rounds in Phase 1 l j length of the j-th frame of Phase 1 k j number of hash functions in the j-th frame of Phase 1 s j random seed used in the j-th frame of Phase 1 U r set of remaining active unexpected tags after Phase 1 N * number of remaining active tags after Phase 1 P 1,j false positive rate in the j-th frame of Phase 1 T 1 time cost of Phase 1 W number of rounds in Phase 2 f w length of the w-th frame of Phase 2 R w number of hash functions in the w-th frame of Phase 2 d w random seed used in the w-th frame of Phase 2 P 2,w false positive rate in the w-th frame of Phase 2 T 2 time taken to execute W rounds in Phase 2 T theoretical execution time q prob. of detect a missing tag in a given slot of Phase 2 Z random variable for slot of the first detection E[T D ] expected detection time of BMTD

Figure 4 . 2 :

 42 Figure 4.2: Example illustrating BMTD

  . 4.2(b), after the first phase, the reader starts to detect missing tags by broadcasting parameters f w and R w . Here we assume W = 1, R w = 2 and f w = 7. By using f w and R w , tag 1 and tag 3 generate a Bloom filter vector, respectively, which is shown in (2) of Fig. 4.2(b). Then they transmit following their individual Bloom filter vector. By sensing the channel, the reader can encode a Bloom filter and use it to check the IDs in E one by one. As shown in (3) of Fig. 4.2(b), since the Bloom filter is constructed based on the responses of tag 1 and tag 3, tag 1 passes the test but tag 2 fails and is regarded as absent. Then the protocol reports a missing event.

Lemma 4 .

 4 3 shows that x * e falls into the range [0, 2x 0 ]. We can thus search [0, 2x 0 ] to find x * e that minimises E[T D ] and then set J and k * j such that ∑ J j=1 k * j = x * e .

Figure 4 . 3 :

 43 Figure 4.3: q pre vs. m.

Fig. 4 .

 4 Fig. 4.3 illustrates the detection probability of SRC with the various number of missing tags under different unexpected tag population sizes. To obtain the figure, we set |E| = 10 3 and = 0.1. It is observed that in the cases that |U| = 0.5 * 10 4 , 1 * 10 4 , 2 * 10 4 , SRC is able to detect at least a missing tag event with probability

Figure 4 . 4 :Fig. 4 .

 444 Figure 4.4:E[T D ] E[T D ] vs. |U| |U| .

Figure 4 . 5 :Figure 4 . 6 :

 4546 Figure 4.5: Actual reliability vs. number of missing tags

Figure 4 . 7 :Figure 4 . 8 :Figure 4 . 9 :Figure 4 . 10 :

 474849410 Figure 4.7: Actual reliability vs. number of unexpected tags

kg ] .

 . The above inequality is readily transformed to the following inequal-

Theorem 5 . 5 .

 55 When f * r in (5.12) and f * gr in (5.21) are power-multiples of two, it holds that T * AB ≥ T * GAB .Proof. Recall Y g in Theorem 5.2 and Z * g in Theorem 5.4 and let x gr k * g n gr , we can rearrange (5.23) as β ir = ∑ G g=1 x gi ∑ G g=1 x gr and φ ir =∑ G g=1 √ x gi x gr , we have: x wr + x gr x wi ),(5.29)

12 Figure 5 . 2 :

 1252 Figure 5.2: Performance comparison in symmetric scenario

12 Figure 5 . 3 :

 1253 Figure 5.3: Performance comparison in asymmetric scenario

  . of missing tags of each group Execution time (slot) B-detect AB-detect GAB-detect (a) Case 1: R = 6, G = 6 . of missing tags of each group Execution time (slot) B-detect AB-detect GAB-detect (b) Case 2: R = 12, G = 6 . of missing tags of each group Execution time (slot) B-detect AB-detect GAB-detect (d) Case 4: R = 12, G = 12

Figure 5 . 4 :

 54 Figure 5.4: Performance comparison in more asymmetric scenario

  This thesis has been dedicated to addressing the fundamental problem of tag counting and monitoring in largescale RFID systems, at both the theoretical modeling and analysis and the algorithm design and optimisation levels, with Chapter 2 focusing on the stability analysis of FSA protocol, Chapter 3 proposing a tag population estimation framework in dynamic RFID systems, Chapter 4 addressing miss tag event detection problem in the presence of unexpected tags, and Chapter 5 devising a suit of algorithms for multiple-group multipleregion RFID systems to detect missing tag event. More specifically, Chapter 2 presented complete and accurate characterisation of FSA behavior, which provides theoretical guidelines on the design of stable FSA-based protocols in other practical applications such as vehicular networks and M2M networks. In Chapter 3, we tackled the dynamic tag population estimation problem and showed a theoretical method of analyzing its estimation error and convergence rate. Furthermore, we illustrated the key to algorithm design and parameter configuration in missing tag event detection problem in Chapter 4 and Chapter 5.

  tions les plus importantes dans les systèmes RFID. Différent des travaux existants dans ce domaine, nous nous concentrons sur deux scénarios inexplorés mais fondamentalement importants, celui en présence d'étiquettes inattendues et celui des systèmes RFID multiple-région multi-groupe. Dans le premier scénario, nous développons un protocole à deux phases à base de filtre de Bloom. Le protocole proposé exploite le filtre Bloom en séquence pour désactiver les étiquettes inattendues d'abord et ensuite teste la composition des étiquettes attendues, limitant ainsi l'interférence des étiquettes inattendues et par conséquent le temps de détection. Pour minimiser le temps de détection tout en garantissant la fiabilité de détection, nous effectuons une analyse théorique et une optimisation sur la configuration des paramètres du protocole. Dans le deuxième scénario, nous formulons et étudions un nouveau problème de détection d'étiquette manquante, survenant dans les systèmes RFID multiple-région multi-groupe, où un lecteur mobile a besoin de détecter s'il y a un événement manquant pour chaque groupe d'étiquettes. L'objectif est de concevoir des protocoles de détection d'étiquette manquante avec le temps d'exécution minimal tout en respectant l'exigence de la fiabilité de détection pour chaque groupe. Nous développons une suite de trois protocoles de détection d'étiquette manquante, chacun diminuant la durée d'exécution par rapport à son prédécesseur en intégrant une version améliorée de la conception du filtre Bloom. En analysant successivement les protocoles développés, nous nous dirigeons progressivement vers un protocole de détection optimal qui fonctionne dans la pratique.

Table 2 .

 2

		1: Main Notations
	Symbols Descriptions
	p	persistence probability
	M	maximum MPR capacity
	Λ	expected arrival rate per slot
	N t	expected arrival rate in frame t
	λ t (n)	prob. of n new arrivals in frame t
	L t	the length of frame t
	X t	random variable: No. of backlogs in frame t
	i	the value of backlogs in frame t, i.e., X t = i
	Y t	random variable: No. of transmitted packet in frame t
	ĥ	the value of packets sent in frame t, i.e., Y t =	ĥ
	Z t	random variable: No. of retransmitted packet in frame t
	h	the value of retransmitted packets in frame t, i.e., Z t = h
	α	the ratio of ĥ to L t
	ξx 0 k 0 ξ ĥk	prob. of having k 0 out of x 0 successful packets in a slot prob. of having k out of ĥ successful packets in frame t
	P is	one-step transition probability
	D

i drift in frame t and its effective operation in practical systems:

• Q1: Under what condition(s) is FSA stable?

Table 3 .

 3 Probability of an idle slot in frame k u k ,Var[u k ] Gaussian random variable and Variance of u k φ k

		1: Main Notations
	z k	System state in frame k: tag population
	y k	Measurement in frame k: idle slot frequency
	ẑk+1|k	Priori prediction of z k+1
	ẑk|k	Posteriori estimate of z k
	P k+1|k	Priori pseudo estimate covariance
	P k|k	Posteriori pseudo estimate covariance
	v k	Measurement residual in frame k
	K k	Kalman gain in frame k
	Q k , R k	Two tunable parameters in frame k
	e k|k-1	Estimation error in frame k
	Φ k	Normalization of v k
	L k	The length of frame k
	Rs k ,h(•)	Random seed in frame k and Hash function
	r k	Persistence probability in frame k
	N k	The number of idle slots in frame k
	p(z k )
		Controllable parameter
	w k	Random variable: variation of tag population
	θ, Υ k	CUSUM threshold and reference value
		Upper bound of initial estimation error
	λ

k , δ k Upper bounds of E[w k ] and E[w 2 k ]

  Definition 3.2 (Boundedness of Random Variable). The stochastic process of the estimation error e k|k-1 is said to be bounded with probability one (w.p.o.), if there exists X > 0 such that (Boundedness in Mean Square). The stochastic process e k|k-1 is said to be exponentially bounded in the mean square with exponent ζ, if there exist real numbers ψ 1 , ψ 2 > 0 and 0

	lim k→∞	k≥1 sup	P{|e k|k-1 | > X} = 0.	(3.11)
	Definition 3.3

.

  The conditions in Lemma 3.2 can be interpreted as follows: To prove the boundedness of e k|k-1 , it is sufficient by constructing a function V k (e k|k-1 ) such that both its drift, i.e,(3.18), and

	V k (e k|k-1 ) k|k-1 e 2	, i.e, (3.16), (3.17), are bounded.

  Lemma 3.5. Under the conditions of Theorem 3.1, if P 0|0 > 0, there exist p

	p	k ≤ P k|k-1 ≤ p k .	(3.50)
	Proof. Recall (3.4) and (3.6), we have P k|k-1 ≥ Q k-1 , and

k

, p k > 0 such that the pseudo-covariance P k|k-1 is bounded for every k ≥ 1, i.e.,

  and Algorithm 2, if there exist time-varying positive real number λ k , σ k > 0 such that Note that the condition E[w k ] ≤ λ k always holds for E[w k ] < 0, we thus focus on the case that E[w k ] ≥ 0. In the proof, the explicit formulas of λ k and σ k are derived. As in the static case, the conditions may be too stringent such that the results still hold even if the conditions are not satisfied, as illustrated in the simulations.The proof of Theorem 3.2 is also based on Lemmas 3.6, 3.7 and 3.8, but due to the introduction of w k into s k , we need another two auxiliary lemmas on E[s k ] and E[s2 k ].

	E[w k ] ≤ λ k ,	(3.71)
	E[w k	2 ] ≤ σ k ,	(3.72)
	then the estimation error e k|k-1 defined by (3.10) is exponentially bounded in mean square and bounded w.p.o..
	Remark.

70) Next, we show the boundedness of the estimation error in Theorem 3.2. Theorem 3.2. Under the conditions of Theorem 1, consider the discrete-time stochastic system given by (3.28) and (3.25) Lemma 3.9. If E[w k ] ≥ 0, then there exists a time-varying real number d k

  .[START_REF] Li | Identifying the missing tags in a large RFID system[END_REF]

	Substituting (3.7), (3.54) and using Lemma 3.8, noticing that E[u k ] = 0, we get
	E	s 2 k P k+1|k	e k|k-1 ≤	1 + φ k φ k P k|k-1	E[w 2 k ] +	2.46 ẑk|k-1
	Lemma 3.10. There exists a time-varying parameter ξ * k > 0 such that E s k 2 P k+1|k	e k|k-1 ≤ ξ * k .
	Proof. By (3.70), we have s 2 k = w 2 k -2K k w k u k + K 2 k u 2 k . Since w k and u k are uncorrelated and e k|k-1 does not
	depend on either w k or u k , we have
		E	s 2 k P k+1|k	e k|k-1 =	E[w 2 k ] P k+1|k	+	K 2 k E[u 2 k ] P k+1|k	.	(3.74)

Table 3 .

 3

					2: Execution time
	Algo.	Variation of tag population 12500 6737 9364 7049 8616 11143 13385 8713 10761
	JREP 3.28	3.28 3.28 3.28 3.28 3.28	3.28	3.28 3.28
	EEKF 1.2	1.2	1.2	0.6	0.6	1.8	0.6	0.6	1.8

  .30) Injecting E[Z] into E[T D ] completes the proof. After deriving E[T D ] as a function of x, we seek the optimum, denoted by x * e , which minimizes E[T D ]. To this end, we first establish an upper-bound of x * e in the following lemma. We write E[T D ] as a function of x. Specifically, let E[T D ] = g(x).To prove the lemma, we show that for any x > 2x 0 it holds that g(x) ≥ g(x 0 ) where x 0

	Lemma 4.3. It holds that x * e ≤ 2t t ln 2 t r |E|q min	.
	Proof. t t ln 2 t r |E|q min	.
	To this end, we first derive the bounds of g(x). Recall (4.27),(4.28), (4.29) and (4.30), we have
		g(x) >	|E|t r x ln 2	,
		g(x) ≤	|E|t r x ln 2	+	t t q min	.

1 .

 1 Set parameters in Phase 2: given |E|, M, α and |U|, compute f * w and y * by (4.17) and (4.18), respectively, and set R w and W such that R w W = y * ; 2. Set parameters in Phase 1: compute x * by Theorem 4.1 if the objective is to minimise the worst-case execution time; compute x * e if the objective is to minimise the expected detection time; then the set of k * j and J is given such that ∑

Table 4 .

 4 2: Actual reliability and detection time of BMTD

	Strategy	10000	Number of unexpected tags 15000 20000 25000	30000
	Worst-1	(1,4108) (1,4441) (1,5013) (1,5453) (1,5510)
	Expected-1	(1,1975) (1,3187) (1,3569) (1,3828) (1,4191)
	Worst-50	(1,1357) (1,1841) (1,2753) (1,2762) (1,2995)
	Expected-50 (1,1353) (1,1618) (1,2272) (1,2472) (1,2815)

Table 5 .

 5

		1: Main Notations
	Symbols Descriptions
	G	Number of groups
	g	Group index and group ID
	R	Number of interrogation regions
	r	Region index
	E	Set of target tags that need to be monitored
	n gr	Number of tags of group g in region r
	m g	Number of missing tags in group g
	M g	Threshold of group g
	P dg	Probability of detecting a missing event of group g
	α g	System requirement on the detection reliability for group g
	f	Length of Bloom filter in B-detect
	k	Number of hash functions in B-detect
	s	Hash function seed
	P f p	False positive rate of Bloom filter in B-detect
	T B	Execution time of B-detect
	f r	Bloom filter vector size in region r in AB-detect
	k g	Number of hash functions for group g in AB/GAB-detect
	P f p,g	False positive rate of Bloom filter for g in AB/GAB-detect
	T AB	Execution time of AB-detect
	f gr	Bloom filter vector length for group g in r in GAB-detect
	T GAB	Execution time of GAB-detect

Phase 2: Virtual Bloom filter construction and missing event detection.

  After interrogating all R regions, the reader combines the Bloom filter vectors B r (1 ≤ r ≤ R) to a virtual Bloom filter B by ORing each bit of them, i.e., B(

5.

3.2 on parameter optimisation. Upon receiving the request, each tag in region r, regardless of the group to which it belongs, selects k slots (h v (ID) mod f )

(1 ≤ v ≤ k)

in the frame of f slots and replies in these slots. The reader then constructs a Bloom filter vector B r with the responses from the tags in each region r as follows. Note there are two types of slots: empty slots and nonempty slots. According to the responses from tags, if slot i (1 ≤ i ≤ f ) is empty, the reader sets B r (i) = 0, otherwise it sets B r (i) = 1.

Phase 2: Virtual Bloom filter construction and missing event detection.

  ≤ • • • ≤ f gR . When receiving the query, each tag compares its group ID with g. If the tag does not belong to the group being queried, it keeps silent and waits for the next query.Otherwise, the tag selects k g positions (h v (ID) mod f gr ) (1 ≤ v ≤ k g ) in the frame of f gr slots and transmits a short response at each of the k g slots. The reader then constructs a Bloom filter for each group g and each region r, denoted by B GAB gr . After interrogating all R regions, the reader combines B GAB gr

  .26) The above optimization problem can be further decomposed to G sub-problem where sub-problem g (1 ≤ g ≤ Solving ∇ f gr ,λ g = 0 yields the following optimum:

	G) is specified as below:
							R
		Minimize:	r=1 ∑	f gr
		Subject to:	R ∑ r=1	n gr f gr	≤ Z g .
	We use the method of Lagrange multiplier to solve each sub-problem g. The Lagrange function can be
	defined as	L( f gr , λ g ) =	R ∑ r=1	f gr + λ g	R ∑ r=1	n gr f gr	-Z g .	(5.27)
		f gr =	√ n gr • ∑ R r=1 Z *	√ n gr
				ln(1-α	1 Mg g	)
					ln 1 2

g

, where Z * g is the maximum of Z g achieved at k * g =

The ergodicity region of a Markov chain in this chapter is referred to as stability region.

The outputs of the hash function have a uniform distribution such that the tag can choose any slot within the round with the equal probability.

For two variables X, Y, asymptotic notation X = Θ(Y) implies that there exist positives c 1 , c 2 and x 0 such that for ∀X > x 0 , it follows that c 1 X ≤ Y ≤ c 2 X.

Acknowledgements

mathematically characterised the performance of the proposed framework in terms of estimation accuracy and convergence speed by deriving the closed-form conditions on the design parameters under which our scheme can stabilise around the real population size with bounded relative estimation error that tends to zero within exponential convergence rate.

Motivation and problem statement

According to the statistics presented in [START_REF]National retail security survey[END_REF], inventory shrinkage, a combination of shoplifting, internal theft, administrative and paperwork error, and vendor fraud, resulted in 44 billion dollars in loss for retailers in 2014. Fortunately, RFID technology can be used to reduce the cost by monitoring products for its low cost and non-line-of-sight communication pattern. Obviously, the first step in the application of loss prevention is to determine whether there is any missing tag. Hence, quickly finding out the missing tag event is of practical importance.

The presence of unexpected tags, however, prolongs the detection time and even leads to missing detection.

Here, we present two examples to motivate the presence of unexpected tags in realistic scenarios.

• Example 1. Consider a retail store with expensive goods and a much larger amount of inexpensive goods, and an RFID system is deployed to monitor the goods. Because of the higher value of expensive products, they are expected to be detected more frequently, but the tags of inexpensive goods also response the interrogation of readers, which influences the decision of readers.

• Example 2. Consider a large warehouse rented to multiple companies where the products of the same company may be placed in different zones according to their individual categories, such as child food and adult food, chilled food and ambient food. When detecting the tags identifying products from one company, readers also receive the feedbacks from the tags of other companies.

In both examples, how to effectively reduce the impact of unexpected tags is of critical importance in missing tag detection. In this chapter, we consider a scenario, as depicted in Fig. 4.1, where each product is affixed by an RFID tag. The reader stores the IDs of expected tags. The problem we address is how to detect

Performance Evaluation

The problem addressed in this chapter is to detect the missing expected tags in the presence of a large number of unexpected tags in a time-efficient and reliable way. In this section, we evaluate the performance of the proposed BMTD. It has been shown in [START_REF] Shahzad | Expecting the unexpected: Fast and reliable detection of missing RFID tags in the wild[END_REF] that existing missing detection protocols cannot achieve the required reliability when there are unexpected tags in the RFID systems except the latest RUN [START_REF] Shahzad | Expecting the unexpected: Fast and reliable detection of missing RFID tags in the wild[END_REF]. We thus compare our proposed BMTD to RUN in terms of the actual reliability and the detection time. Note that the detection time can be interpreted as the time taken to either detect the fist missing tag event if a missing tag is found or complete the execution if no missing tag is found.

The simulation parameters are set with reference to [START_REF] Luo | Missing-tag detection and energy-time tradeoff in large-scale RFID systems with unreliable channels[END_REF] and [START_REF] Shahzad | Expecting the unexpected: Fast and reliable detection of missing RFID tags in the wild[END_REF]. Specifically, since both transmission rates from the tags to the reader and the reader to the tags depend on physical implementation and interrogation environment, we make the same assumption as in [START_REF] Luo | Missing-tag detection and energy-time tradeoff in large-scale RFID systems with unreliable channels[END_REF] that t r = t t . Moreover, because RUN is the baseline protocol, we use the similar simulation scenarios and the same performance metrics as in [START_REF] Shahzad | Expecting the unexpected: Fast and reliable detection of missing RFID tags in the wild[END_REF] where the time needed to detect a missing tag event is shown in terms of the number of slots. To that end, we, without loss of generality, assume t r = t t = 1 in (4.26) in the simulation. Besides, we compute the optimal parameter values for RUN by following its specifications.

In the simulation, we use SRC [START_REF] Chen | Understanding RFID counting protocols[END_REF] armed with missing tag detection function in this chapter to estimate the unexpected tag population size with the confidence rang = 0.1. And all presented results are obtained by taking the average value of 100 independent trials under the same simulation setting.

We start by evaluating the performance of the BMTD by optimizing the worst-case execution time and the expected detection time.

Comparison between two strategies of BMTD

In this subsection, we compare the performance of two strategies of the BMTD which are abbreviated to Worst-M and Expected-M here, respectively. We set |E| = 1000, m = 100, α = 0.9, |U| = 10000 : 5000 : 30000, M = 1 and 50. Table 4.2 lists the results where the first and second elements in the two-tuple (•, •) denote the actual reliability and detection time, respectively. It can be seen that Expected-M costs less time than Worst-M to achieve the same reliability which is greater than the system requirement on the detection reliability, especially when M is small. Specifically, compared with Worst-1, Expected-1 reduces the detection time by up to 51.92% when |U| = 10000. This is because x * = 5 is too large for Phase 1 by optimizing the worst-case execution time, which wastes time. In contrast, minimizing the expected detection time relieves the influence of unexpected tag population size on the time of Phase 2 and thus outputs a smaller x * e = 2. In the rest of our simulation, we configure the parameters of the BMTD to minimise the expected detection time. tags in order to reduce their interference to the detection. In the second phase, we further used Bloom filter to test the membership of the expected tags to detect missing tags. We also showed how to configure the protocol parameters so as to optimize the detection time with the required reliability. Furthermore, we conducted extensive simulation experiments to evaluate the performance of the proposed protocol and the results demonstrate the effectiveness and efficiency of the propose protocol in comparison with the state-of-the-art solution.

Chapter 5

On Missing Tag Detection in Multiple-group

Multiple-region RFID Systems

Introduction

We investigate a different version of missing tag detection problem in this chapter motivated by the increasing application of mobile reader [START_REF] Liu | Towards adaptive continuous scanning in large-scale rfid systems[END_REF][START_REF] Xie | Exploring the gap between ideal and reality: An experimental study on continuous scanning with mobile reader in RFID systems[END_REF] and the following practical settings.

• Multiple groups of tags. Tags are usually attached to objects belonging to different groups: e.g., different brands of the goods with the high-end brands order-of-magnitude more valuable than their low-end peers.

Therefore, the missing tag events are characterized by asymmetrical threshold and reliability requirement across groups.

• Multiple interrogation regions. Tags may be unevenly located in multiple interrogation regions: e.g., tags may be located in several rooms or different corners or regions of a large warehouse. Hence, a reader may need to move several times to cover all monitored tags and complete the missing tag detection process.

The problem we consider is to devise missing tag detection protocol with minimum execution time while guaranteeing the detection reliability requirement for each group of tags in multiple-region scenario. In the considered multiple-group multiple-region scenario, all existing missing tag detection protocols cannot work effectively due to the following two reasons. First, existing approaches as reviewed in Section 4.2 of Chapter 4 require the full coverage of tags when executing the detection algorithms, which clearly does not hold in the considered multiple-region scenario. Secondly, existing work does not take into account the heterogeneity among groups and thus either cannot meet the individual reliability requirement, or suffers extremely long detection delay.

To solve this challenging problem, we deliver a comprehensive analysis on the missing tag detection problem in the above multiple-group multiple-region environment and investigate how to devise optimum missing tag detection algorithms. Note that when there are only one group and all tags are with one interrogation region, our problem degenerates to the classical missing tag detection problem studied in the literature.

To design missing tag detection algorithms in the multiple-region multiple-group case, we leverage a powerful technique called Bloom filter which is a space-efficient probabilistic data structure for representing a set and supporting set membership queries [START_REF] Bloom | Space/time trade-offs in hash coding with allowable errors[END_REF] to detect a missing event. Specifically, we develop a suite of three missing tag detection protocols, each decreasing the execution time compared to its predecessor by incorporating an improved version of the Bloom filter design and parameter tuning. By sequentially analysing the developed protocols, we gradually iron out an optimum detection protocol that works in practice.

System Model and Problem Formulation

System Model

We consider a grouped RFID system composed of a mobile reader and G groups of tags distributed in R (R ≥ 1) interrogation regions (e.g., R rooms), concisely referred to as regions. In case where a tag may be physically located in two regions, i.e., regions may overlap one with another, the tag only responses to reader queries regarding to the first region when it is interrogated. In this sense, we can treat the regions as non-overlapping ones.

We use E to denote the set of the tags which are expected to be present and we denote its cardinality (i.e., the number of expected tags) by |E|. The reader knows the IDs of all tags in E but does not know the set of tags in each region. For presentation conciseness, we set the ID of group g (1 ≤ g ≤ G) to its index g. We assume every tag knows its group ID through a grouping protocol, e.g. [START_REF] Liu | Fast rfid grouping protocols[END_REF]. We also assume the reader knows the approximate number of tags of each group g actually present in each region r (1 ≤ r ≤ R), denoted by n gr .

The estimation of n gr can be achieved by the reader by deactivating all tags not belonging to group g (using the ID of group g) and then using any state-of-the-art tag population estimation algorithm. Table 5.1 summaries main notations used in this chapter.

Problem Formulation

We are interested in detecting missing tag event for each group g. Let m g denote the number of missing tags in group g which is of course not known by the reader. Let M g denote the threshold of group g. A missing event of group g denotes the event where there are at least M g tags of group g missing in the system. Let P dg denote the probability that the reader can detect a missing event of group g, we formulate the optimum missing tag detection problem as follows.

Definition 5.1 (Optimum missing tag detection problem). The optimum missing tag detection problem is to devise an algorithm of minimum execution time which can detect a missing event for each group g with probability P dg ≥ α g if m g ≥ M g , where α g is the requirement on the detection reliability for group g. When there is only one group in the system, the problem degenerates to the classical missing event detection problem.

It can be easily checked that f achieves its minimum as (5.4) at k * . The theorem is thus proved.

Remark. Given the practical meaning of k * and f * , both of them should been further rounded to the smallest integers not smaller than themselves.

The Adaptive Approach

In B-detect, the reader uses the same parameters in each region, particularly the length of the Bloom filter, which may not be optimum in the case when the tags are not evenly distributed across interrogation regions.

Motivated by this observation, we develop another missing tag detection protocol, named AB-detect, which enables the reader to use different parameters based on the number of tags in the region the reader queries.

Protocol Description

Phase 1: Query and feedback collection. The reader performs a query in each region r with the parameter

, s) where f r is the length of the Bloom filter vector used in region r, k g is the number of hash functions used by tags in group g, s is the hash seed which is identical for all groups and regions. There are two differences compared to the baseline approach. First, f r may be different across different regions but identical across groups; Second, k g may be different across different groups but identical across regions. We require f r to be a power-multiple of two, i.e., f r = 2 b r , (b r ∈ N). As in B-detect, the reader constructs an f r -bit Bloom filter vector B r with the responses from the tags in each region r. Without loss of generality, we assume

as illustrated in Fig. 5.1. The reader then performs membership test. For each tag in group g, the reader maps its ID into k g bits at positions (h v (ID)

If all the corresponding bits in B AB are 1, then the tag is regarded as present. Otherwise, the tag is considered to be missing. The reader reports a missing event for group g if the number of missing tags in the group g is at least M g and no missing event otherwise. The following lemma proves that there is no false negative in AB-detect.

When the detection requirement is not stringent, e.g., there are a large number of missing tags and the reader only needs to detect a small number of them so as to report a missing event, estimating the number of tags may be used. However, when the detection requirement is stringent, estimating the number of tags is not efficient as it either requires long execution time or cannot satisfy the detection requirement. To demonstrate this, we have conducted more experiments by comparing our approach with the estimation of tag numbers. Under the same detection reliability requirement, the estimation algorithm spends over 48 -72.6 times as much time as our algorithms. Therefore, in our approach, we perform a coarse estimation on the tag population for two reasons: 1) our algorithms need a coarse estimation of tag population to configure parameters; 2) in case when the detection requirement is not stringent, this phase allows the reader to quickly detect a missing event.

Presence of Unknown/Unexpected Tags

Unknown and unexpected tags can be interpreted as the tags that have not been identified by the reader [START_REF] Liu | Unknown tag identification in large RFID systems: An efficient and complete solution[END_REF],

such as newly arrived products, on which the reader does not have any knowledge. During the interrogation, the unknown tags will respond together with the known tags, which results in the interference to the detection of missing known tags and thus degrades the performance [68] [91].

Fortunately, two of our proposed algorithms, AB-detect and GAB-detect, are resistant to the interference caused by unknown tags. The reason is as follows. The unknown tags have not been identified by the reader, so they do not have their individual group IDs [START_REF] Liu | Fast rfid grouping protocols[END_REF] such that no group ID in the interrogation messages matches with theirs. Therefore, unknown tags stay silent during the whole detection process.

Numerical Results

In this section, we evaluate the performance of the proposed approaches in terms of execution time and investigate tradeoffs under different parameter settings.

Simulation Settings

We conduct the experiment under both symmetric and asymmetric scenarios under different settings of R, G and M g . By symmetric/asymmetric, we mean that tag population size in each region r is identical/different.

Moreover, we set the same M g for all group g but vary α g for each group. Moreover, we use the symmetric transmission rate as in [START_REF] Chen | simples à formuler et intuitivement compréhensible, tandis que tous les deux présentent des défis importants à la fois fondamentaux et pratiques, et exigent des efforts non négligeables à résoudre. Plus précisément, nous abordons les problèmes suivants allant de la modélisation et de l'analyse théorique[END_REF] [72] in the numerical analysis and set the transmission time for one bit to be one slot, i.e., δ = 1. The length of group ID is set to log 2 G bits as in [START_REF] Liu | Fast rfid grouping protocols[END_REF]. We simulate the optimum parameters settings derived in (5.4) (5.5) for B-detect, (5.12) (5.13) for AB-detect, and (5.21) (5.22) for GAB-detect.

For a comprehensive evaluation, we simulate four cases with different combination of (R, G) in both the symmetric and asymmetric scenarios: case 1: (6, 6), case 2: [START_REF] Liu | Efficient unknown tag identification protocols in large-scale rfid systems[END_REF][START_REF] Noh | Anc-aloha: Analog network coding ALOHA for satellite networks[END_REF], case 3: [START_REF] Noh | Anc-aloha: Analog network coding ALOHA for satellite networks[END_REF][START_REF] Liu | Efficient unknown tag identification protocols in large-scale rfid systems[END_REF], and case 4: [START_REF] Liu | Efficient unknown tag identification protocols in large-scale rfid systems[END_REF][START_REF] Liu | Efficient unknown tag identification protocols in large-scale rfid systems[END_REF]. The required detection reliability for group g (1 ≤ g ≤ G) is set to α g = 0.749 + 0.05(g -1), i.e., 0.749≤α g ≤0.999 in group g. Moreover, we also investigate the impact of estimation error on the performance.

From the results listed in Table 5.2, we can observe that GAB-detect significantly outperforms AB-detect and B-detect when M g is different for each group. Besides, the execution time increases by less than 11% in the worst case when varies from 0 to 10%. While on average, B-detect and GAB-detect and AB-detect spend 9% and 6% and 2.6% more time, respectively. Therefore, it is fair to allow = 10%.

Conclusion

In this chapter, we formulate a missing tag detection problem arising in multiple-group multiple-region R-FID systems, where a mobile reader needs to detect whether there is any missing event for each group of In this thesis, we present a systematic research on a number of research problems related to tag counting and monitoring, one of the most fundamental component in RFID systems, particularly when the system scales. These problems are simple to state and intuitively understandable, while of both fundamental and practical importance, and require nontrivial efforts to solve. Specifically, we address the following problems ranging from theoretical modeling and analysis, to practical algorithm design and optimization. To lay the theoretical foundations for the algorithm design and optimization, we start by studying the stability of frame slotted Aloha. We model system backlog as a Markov chain. The main objective is translated to analyze the ergodicity of the Markov chain. We then establish a framework of stable and accurate tag population estimation schemes based on Kalman filter for both static and dynamic RFID systems. We further proceed to addressing the problem of missing tag detection, one of the most important RFID applications. Different from existing works in this field, we focus on two unexplored while fundamentally important scenarios, missing tag detection in the presence of unexpected tags and in multiple-group multiple-region RFID systems. In the first scenario, we develop a Bloom filter-based protocol. In the second scenario, we develop three protocols by incorporating an improved version of the Bloom filter design.