Wenceslao Palma

Reza Akbarinia
email: rakbarin@cs.uwaterloo.ca

Esther Pacitti

Patrick Valduriez
email: patrick.valduriez@inria.fr

Efficient Processing of Continuous Join Queries using Distributed Hash Tables

This paper addresses the problem of computing approximate answers to continuous join queries. We present a new method, called DHTJoin, which combines hash-based placement of tuples in a Distributed Hash Table (DHT) and dissemination of queries using a gossip style protocol. We provide a performance evaluation of DHTJoin which shows that DHTJoin can achieve significant performance gains in terms of network traffic.

Introduction

Recent years have witnessed major research interest in data stream management systems. A data stream is a continuous and unbounded sequence of data items. There are many applications that generate streams of data including financial applications [START_REF] Chen | NiagaraCQ: A Scalable Continuous Query System for Internet Databases[END_REF], network monitoring [START_REF] Sullivan | Tribeca: A System for Managing Large Databases of Network Traffic[END_REF], telecommunication data management [START_REF] Chandrasekaran | TelegraphCQ: Continuous Dataflow Processing for an Uncertain World[END_REF], sensor networks [START_REF] Bonnet | Towards Sensor Database Systems[END_REF], etc. Processing a query over a data stream involves running the query continuously over the data stream and generating a new answer each time a new data item arrives. Due to the unbounded nature of data streams, it is not possible to store the data entirely in a bounded memory. This makes difficult the processing of queries that need to compare each new arriving data with past ones. We are interested in systems which have limited main memory but that can tolerate an approximate query result which has a maximum subset of the result. An example of such queries is join queries which are very important for many applications. As an example, consider a network monitoring application that needs to issue a join query over traffic traces from various links, in order to monitor the total traffic that is common among three links L1, L2 and L3 over the last 10 minutes. Each link (stream) contains tuples each one with a packet identifier pid and the packet size. This query can be posed using a declarative language such as CQL [START_REF] Arasu | An Abstract Semantics and Concrete Language for Continuous Queries over Streams and Relations[END_REF], which is a relational query language for data streams, as follows:

Select sum (L1.size) From L1[range 10 min], L2[range 10 min], L3[range 10 min] Where L1.pid=L2.pid and L2.pid=L3.pid A common solution to the problem of processing join queries over data streams is to execute the query over a sliding window [START_REF] Golab | Processing Sliding Windows Multi-Joins in Continuous Queries over Data Streams[END_REF] that maintains a restricted number of recent data items. This allows queries to be executed in a finite memory and in an incremental manner by generating new answers when a new data item arrives.

In this paper, we address the problem of computing approximate answers to windowed stream joins over data streams. Our solution involves a scalable distributed sliding window that takes advantage of the free computing power of DHT networks and can be equivalent to thousands of centralized sliding windows. Then, we propose a method, called DHTJoin, which deals with efficient execution of join queries over all data items which are stored in the distributed sliding window. DHTJoin combines hash-based placement of tuples in the DHT and dissemination of queries using a gossip style protocol. We evaluated the performance of DHTJoin through simulation. The results show the effectiveness of our solution compared with previous work.

The rest of this paper is organized as follows. In Section 2, we introduce our system model and define the problem. In Section 3 we describe DHTJoin. Section 4, describes a performance evaluation of our solution through simulation using Simjava. In Section 5, we discuss related work. Finally, Section 6 concludes.

System Model and Problem Definition

In this section we introduce a general system model for processing data streams over DHTs, with a DHT model, a stream processing model and a gossip dissemination system. Then, we state the problem.

DHT Model

In our system, the nodes of the overlay network are organized using a DHT protocol. While there are significant implementation differences between DHTs [START_REF] Ratnasamy | A Scalable Content-Addressable Network[END_REF] [START_REF] Stoica | Chord: A Scalable Peer-to-Peer Lookup Service for Internet Applications[END_REF], they all map a given key k onto a peer p using a hash function and can lookup p efficiently, usually in O(log n) routing hops where n is the number of peers. DHTs typically provide two basic operations [START_REF] Stoica | Chord: A Scalable Peer-to-Peer Lookup Service for Internet Applications[END_REF]: put(k, data) stores a key k and its associated data in the DHT using some hash function; get(k) retrieves the data associated with k in the DHT. Tuples and continuous queries are originated at any node of the network. Nodes insert data in the form of relational tuples and each query q is represented in SQL. Tuples and queries are timestamped to represent the time that are inserted in the network by some node. Additionally, each query is associated with a unique key used to identify it in query grouping and to relate it to the node that submitted it.

Stream Processing Model

A data stream S i is a sequence of tuples ordered by an increasing timestamp where i [1..m] and m≥2 denotes the number of input streams. At each time unit, a number of tuples of average size l i arrives to stream S i . We use  i to denote the average arrival rate of a stream S i in terms of tuples per second.

Many applications are interested in making decisions over recently observed tuples of the streams. This is why we maintain each tuple only for a limited time. This leads to a sliding window W i over S i that is defined as follows. Let T(W i ) denotes the size of W i in terms of seconds, i.e. the maximum time that a tuple is maintained in W i . Let TS(s) be a function that denotes the arrival time of a tuple sS, and t be current time. Then W i is defined as

W i = {s| sS i  t-TS(s)  T(W i )}.
Tuples continuously arrive at each instant and expire after T(W i ) time steps (time units). Thus, the tuples under consideration change over time as new tuples get added and old tuples get deleted. In practice, when arrival rates are high and the memory dedicated to the sliding window is limited, it becomes full rapidly and many tuples must be dropped before they naturally expire. In this case, we need to decide whether to admit or discard the arriving tuples and if admitted, which of the existing tuples to discard. This kind of decision is made using a load-shedding strategy [START_REF] Kang | Evaluating windows joins over unbounded streams[END_REF][21] which yields that only a fraction of the complete result will be produced.

Gossip Dissemination System

The basic idea behind gossiping [15][26] is to have all nodes collaborating to disseminate information into the entire network using a partial view stored in a local list of size L. To this end, when a node wishes to disseminate a message or receive it for the first time, it picks k nodes from its local list and sends them the message. Initializing and maintaining the local list at each node in face of dynamic changes in the network is done by a membership protocol. The number of gossip targets, k, is a typical configuration parameter called fanout. The sum of the k links between nodes specifies an overlay network on top of the existing network topology.

If a node receives a message for the second time, it simply discards the message. A relatively inexpensive optimization is to never forward a message back to the node it was just received from. Due to their inherent redundancy, gossip dissemination systems are able to mask network omissions and also node failures. In order to evaluate the performance of gossip dissemination systems, it is essential to define a set of metrics:

 Hit Ratio: the number of nodes that receive a sent message. Ideally, a gossip dissemination system should always achieve a hit ratio of 100%, that is, to reach every node in the network.  Dissemination Speed: the time a message requires to reach every node in the network. It depends of network latency and the number of hops a message takes to reach the last node. In our evaluation we focus on the latter factor.  Message Overhead: the number of redundant times that a message is forwarded during its dissemination. An effective dissemination system assures a high Hit Ratio, a fast dissemination speed, and has a low message overhead.

Problem Definition

In this paper, we deal with the problem of processing join queries over data streams. A data stream is a sequence of tuples ordered by monotonically increasing timestamps. The timestamps are generated using the KTS service [START_REF] Akbarinia | Data currency in replicated DHTs[END_REF] which is a distributed service that deals with generating monotonically increasing timestamps in DHTs. Each tuple and query have a timestamp that may either be implicit, i.e. generated by the system at arrival time, or explicit, i.e. inserted by the source at creation time. However, we do not include the timestamp attribute in the schema of the stream. Formally, the problem can be defined as follows. Let ={S 1 ,S 2 ,……, S m } be a set of relational data streams, and Q={Q 1 ,Q 2 ,....,Q n } be a set of join queries specified on these data streams. Our goal is to provide an efficient method to execute Q over  in terms of network traffic.

DHTJoin Method

Using structured overlay P2P networks in data stream processing environments is a novel research topic [START_REF] Idreos | Continuous Multi-Way Joins over Distributed Hash Tables[END_REF]. We can take advantage of it strengths to support queries over data streams [START_REF] Babcock | Models and Issues in Data Streams[END_REF]. The main issues for processing continuous queries in DHTs are the following: how to route data and queries to peers in an efficient way; how to provide a data storage mechanism for storing relational data; and how to provide a good approximate answer to join queries.

We describe our solution for continuous join query processing using DHTs. Our method has two steps: indexing of tuples and dissemination of queries. A tuple inserted by a node is indexed, i.e., stored at another node using DHT primitives only if there is a query that requires it. To assure the knowledge of that query to the entire network, it is disseminated using a gossip protocol.

Indexing Tuples

Existing DHT implementations such as Chord [START_REF] Stoica | Chord: A Scalable Peer-to-Peer Lookup Service for Internet Applications[END_REF] are robust: they use simple algorithms with provable properties even under frequent node failures, departures and re-joins. To answer a query, we consider different kinds of peers. The first kind is Stream Reception Peers (SRP) which are responsible for disseminating queries and indexing tuples to the second kind of peers named Stream Query Peers (SQP). SQPs are responsible for executing queries over the arriving tuples using their local sliding windows, and sending the results to the third kind of node(s) named User Query Peer (UQP). To support dissemination of queries each node is a gossip node and for indexing of tuples each node is a DHT peer. Note that the distribution between SRP, SQP and UQP is functional and the same peer can have all of these functionalities.

The tuples arrive at any SRP and are indexed onto SQPs if there exists a query that requires it. A join processing algorithm is executed at each SQP that receives the tuples sent by an SRP. In our solution, the tuples are indexed in SQPs. Each arriving tuple is entirely indexed in an SQP using a join attribute value as storage key, usually in O(log n) routing hops where n is the number of peers.

Let us describe our indexing method for three streams S 1 , S 2 and S 3 . However, our method is not limited to three streems and works on any number of streams. Let A be the set of attributes of S 1 , S 2 and S 3 . Let s 1 , s 2 , and s 3 be tuples belonging to S 1 , S 2 and S 3 respectively and val(s a) be a function that returns the value of an attribute a A in tuple s. Let h be a uniform hash function that hashes its inputs into a DHT key, i.e. a number which can be mapped by the DHT onto a peer. For indexing a tuple s 1 that arrives at a SRP, each tuple obtains an index key computed as key= h(S 1 ,val(s 1 ,a)) and is then indexed in a SQP using put(key,s 1 ) only if there exists a query that requires this tuple indexed by using the value of attribute a. For s 2 and s 3 tuples, we proceed in the same way. All tuples of S 1 , S 2 and S 3 having the same value in attribute a get indexed at the same SQP. Join operations are performed on the tuples stored within sliding windows. Let W i denote a sliding window on stream S i . We use time-based sliding windows where T(W i ) is the size of the window in time units. At time t, a tuple belongs to W i if it has arrived in the time interval [t-T(W i ),t].

To illustrate, let us consider an equijoin query over three streams S 1 , S 2 and S 3 . Once a S 1 -tuple arrives onto a SRP we verify if a query requires this tuple indexed by some attribute, if so, we index the tuple and execute a Join processing algorithm at each SQP that receives SRP tuples. Processing a join query over timestamp based sliding windows with three input relational streams S 1 , S 2 and S 3 is done as follows.

Upon each arrival of a new tuple from stream S 1 to W 1 , expired tuples in W 2 and W 3 are invalidated, then for each W 2 -tuple we evaluate the join of W 3 with the arriving tuple and probe the result set. A load shedding procedure is executed over W 1 's buffer if there is not enough memory to insert the tuple. The UQP receives the tuples from SQPs. At any time, a SQP could contain a collection of tuples representing a portion of the join query result. The join query result is inserted into a local queue and sent to the peer that submitted the query. The output of the join sent to that peer (an UQP node) consists of all pairs of tuples s 1  S 1 , s 2  S 2 and s 3  S 3 such that s 1 .a= s 2 .a=s 3 .a and, a time t, both s 1  W 1 , s 2  W 2 and s 3  W 3 .

Disseminating Queries

A query q can originate at any of the nodes and is disseminated to the whole network. This dissemination allows a global knowledge of queries and to index tuples using an attribute value only if there exists a query that requires it. However, disseminating q can be difficult if nodes have only limited knowledge of the members of the network. This knowledge is stored in a local list of size L and we must assure a Hit Ratio of 100%. We assume a network of N nodes. Intuitively, the requirement for assuring a Hit Ratio of 100% is that the components of the local list form a strongly connected directed graph over all nodes. For example, if we use the Chord protocol we can build a local list of size 1 with an immediate successor pointer from the finger table and disseminate a query in N hops without redundant messages. However, the dissemination speed is very slow. To accelerate the dissemination process we can consider a local list composed by an immediate successor pointer from the finger table and the predecessor pointer, used to walk counter clockwise around the Chord ring, maintained by the Chord protocol. In this case, the dissemination speed is reduced to N/2 hops with 2 redundant messages. Considering that a tuple can be indexed in log(N) hops we must increase the dissemination speed as high as possible.

To this end, we can increase the fanout by exploiting the entire finger table and increase the dissemination speed to log(N) hops, i.e., we disseminate a query as fast as a tuple index but that implies redundant messages. Redundancy provides tolerance to dynamic environments but excessive redundancy leads to excessive network traffic.

In addition, with predecessor and successor pointers of the Chord protocol we can form a strongly bidirectional connected directed graph that assures a 100% Hit Ratio and with the entire finger table, we can increase the dissemination speed as fast as a tuple index.

Optimization

Achieving a dissemination speed of queries equal to log(N) hops (taking advantage of the Chord finger table) yields an increase in the number of forwarded messages. We propose an optimization to decrease considerably the number of redundant messages without decreasing the dissemination speed and guaranteeing a 100% Hit Ratio. Our solution is as follows. Each time a node originates a query, it builds a gossip message that consists of the query, its own identifier and its gossip targets: the successor and predecessor pointers and its finger table. Upon reception of a new gossip message, a node builds its gossip targets considering the intersection of the nodes received with the arriving message, its own identifier and its own gossip targets. Thus, we can decrease the number of redundant messages to 1/2log(N) (see our performance evaluation), while maintaining the same dissemination speed and assuring a 100% of Hit Ratio.

Performance Evaluation

To test our DHTJoin method, we implemented a simulator using SimJava [START_REF] Simjava | [END_REF] running on a Linux PC with a 3.4 Ghz Pentium Processor and 512 megabytes of RAM. Our simulator is based on Chord which is a simple and efficient DHT. To simulate a peer, we use a SimJava entity that performs all tasks that must be done by a peer in the DHT, in the dissemination system and in the stream window-based join processing.

We generate arbitrary input data streams consisting of synthetic asynchronous data items with no tuple-level semantics. We have a schema of 4 relations, each one with 10 attributes. In order to create a new tuple we choose a relation and assign values to all its attributes using a Zipf distribution with a default parameter of 0.9. The max value of the domain of the join attribute is fixed to 1000. Tuples on streams are generated at a constant rate of 30 tuples per second. Queries are generated with a mean arrival rate of 0.02, i.e., a query arrives every 50 seconds on average. The network size is set to 1024 nodes. In all experiments, we use time-based sliding windows of 50 seconds. For each of our tests, we run the simulator for 300 seconds. In order to assess our approach, we compare the network traffic of DHTJoin against a complete implementation of RJoin [START_REF] Idreos | Continuous Multi-Way Joins over Distributed Hash Tables[END_REF] which is the most relevant related work (see Section 5). RJoin uses incremental evaluation based on tuple indexing and query rewriting over distributed hash tables. In RJoin a new tuple is indexed twice for each attribute it has; wrt the attribute name and wrt the attribute value. A query is indexed waiting for matching tuples. Each arriving tuple that is a match causes the query to be rewritten and reindexed at a different node. The network traffic is the total number of messages needed to index tuples and disseminate a query in DHTJoin or to index tuples and perform query rewriting in RJoin.

In order to show the effectiveness of the dissemination of queries, we use the Hit Ratio, Dissemination Speed and Message Overhead metrics. The fanout parameter is set to 2, m/2, m and m OPT , with m equal to log(N). The fanout value m OPT includes the optimization process presented in Section 3.

In the rest of this section, we evaluate network traffic and the effectiveness of gossip dissemination.

Network Traffic

In this section, we investigate the effect of tuples' arrival rate, query's arrival rate and number of joins on the network traffic. The network traffic of RJoin and DHTJoin grows as the tuples' arrival rate grows. In RJoin, as more tuples arrive, the number of messages related to the indexing of tuples and query rewriting increases (see Figure 1). As expected, DHTJoin generates significantly less messages because a high tuples' arrival rate does not mean more indexing of tuples. The reason is that before indexing a tuple, DHTJoin checks for the existence of a query that requires it. In Figure 2, we show that, as more queries arrive, RJoin generates more query rewriting messages. However, DHTJoin generates more messages only if new queries related to new different values used in the tuples arrive in the system. Figure 3 shows that more join require more network traffic. RJoin generates more query rewriting when there are more joins in the queries. However, in DHTJoin the network traffic increases only if the arriving queries require of attributes not present in the other queries. The reason is that with the dissemination of queries, DHTJoin can avoid the unnecessary indexing of tuples that are not required by the queries In summary, due to the integration of query dissemination and hash-based placement of tuples our approach avoids the excessive traffic generated by RJoin which is due to its method of indexing tuples.

Gossip Dissemination

To evaluate the effectiveness of gossip dissemination of queries, we post a query originated randomly at any node. In this experiment, we show the dissemination speed and look at the evolution of dissemination in terms of the number of nodes that

have not yet been contacted as a function of the hops taken. Figure 4 shows that our method yieds a complete dissemination, i.e., a Hit Ratio of 100% in log(N) hops. This is because the nodes of the local list form a strongly connected directed graph over all nodes. Our proposed optimization to reduce message overhead produces no changes in the speed of the dissemination. If a query is disseminated as fast as possible we avoid excessive store time of tuples waiting for a query and we can generate join tuples as early as possible. Thus, we state that a query must be disseminated as fast as a tuple. In Figure 5, we confirm this and we can see that the higher the fanout the dissemination speed increases to reach the same speed as a tuple index. However, if we choose a higher fanout, we incur message overhead which overload the network. Conversely, we can choose a lower fanout to obtain a minimal overhead, as shown in Figure 6, but this incurs slower dissemination, as shown in Figure 5. To reduce message overhead we ran experiments to include the optimization proposed in Section 3. In Figure 5, we show that with a fanout m OPT , we obtain a fast dissemination and at the same time, as shown in Figure 6, we reduce the message overhead to 1/2log(N) wrt a fanout m.

Related Work

A DHT can serve as the hash table that underlies many parallel hash-based join algorithms. However, our approach provides Internet-level scalability. Our work is related to many studies in the field of centralized and distributed continuous query processing [START_REF] Huebsch | Queying the Internet with PIER[END_REF][10] [START_REF] Tatbul | Window-Aware Load Shedding for Aggregations Queries over Data Streams[END_REF][6] [START_REF] Liu | A dynamically Adaptive Distributed System for Processing Complex Continuous Queries[END_REF]. In PIER [START_REF] Huebsch | Queying the Internet with PIER[END_REF], a query processor is used on top of a DHT to process one-time join queries. Recent work on PIER has been developed to process only continuous aggregation queries. PeerCQ [START_REF] Gedik | PeerCQ: A Descentralized a Self-Configuring Peer-to-Peer Information Monitoring System[END_REF] was developed to process continuous queries on top of a DHT, However, PeerCQ does not consider SQL queries and the data is not stored in the DHT. Borealis [START_REF] Tatbul | Window-Aware Load Shedding for Aggregations Queries over Data Streams[END_REF], TelegraphCQ [START_REF] Chandrasekaran | TelegraphCQ: Continuous Dataflow Processing for an Uncertain World[END_REF] and DCAPE [START_REF] Liu | A dynamically Adaptive Distributed System for Processing Complex Continuous Queries[END_REF] have been developed to process continuous queries in a cluster setting and many of their techniques for load-shedding and load balancing are orthogonal to our work. The most relevant previous work regarding the utilization of a structured overlay P2P network is [START_REF] Idreos | Continuous Multi-Way Joins over Distributed Hash Tables[END_REF] which proposes RJoin, an algorithm that uses incremental evaluation. This incremental evaluation is based on tuple indexing and query rewriting over distributed hash tables. A major difference in our work differs is that our tuple index mechanism indexes tuples only if there exists a query that requires it. Large-scale information dissemination based on gossip protocols are essential for applications such as replicated database maintenance [START_REF] Demers | Epidemics Algorithms for Replicated DB Maintenance[END_REF], publish/subscribe systems [START_REF] Eugster | Lightweight Probabilistic Broadcast[END_REF] and distributed failure detection [START_REF] Van Renesse | A Gossip-Style Failure Detection Service[END_REF]. These approaches do not consider a structured overlay setting. The probabilistic dissemination algorithm named Randcast proposed in [START_REF] Kermarrec | Probabilistic Realiable Dissemination in Large-Scale Systems[END_REF] spreads messages very fast but fails to reach every node in the network. In [START_REF] Patel | Achieveing Predictable Gossip Dissemination by Leveraging Social Network Principles[END_REF] the authors propose a technique that leverages simple social network principles enabling nodes to select gossip targets intelligently. This work shows experimentally that gossip dissemination on the Chord overlay fails to deliver the messages to every node in the network. Our work instead shows that it is possible to reach every node in the network forming a strongly bidirectional connected directed graph composed by predecessor and successor pointers. In [START_REF] Lin | Gossip versus Deterministic Flooding: Low Message Overhead and High-Reliability for Broadcasting on Small Networks[END_REF] the authors assure a good tradeoff between message overhead and reliability guarantee using a connection graph called a Harary graph. However, its principal drawback is the maintenance of such graph that requires global knowledge of membership. In our work, the structure that supports the membership protocol is supported by the structured overlay and does not require global knowledge of membership for its maintenance. The most relevant previous work regarding dissemination is the hybrid dissemination proposed in [START_REF] Voulgaris | Hybrid Dissemination: Adding Determinism to Probabilistic Multicasting in Large-Scale P2P Systems[END_REF]. This works uses deterministic links to assure a complete dissemination. However, our work differs in the presence of a structured overlay and a membership protocol based on its routing table.

Conclusion

In this paper, we proposed a new method, called DHTJoin, for processing continuous join queries using DHTs. DHTJoin combines hash-based placement of tuples and dissemination of queries using a gossip style protocol. Our performance evaluation shows that DHTJoin obtains significant performance gains due to our schema of global knowledge of queries based on a gossip protocol. This schema has a low message overhead and avoids the excessive traffic produced by the tuple index method and the query rewriting of RJoin. Our results demonstrate that the total number of messages transmitted by DHTJoin is always fewer that RJoin wrt tuple arrival rate, query arrival rate and number of joins.

As future work, we plan to address the problem of efficient execution of top-k join queries over data streams using DHTs taking advantage of the best position algorithms [START_REF] Akbarinia | Best Position Algorithms for Top-k Queries[END_REF] which can be used in many distributed and P2P systems for efficient processing of top-k queries.

Figure 3

 3 Figure 3 Effect of Number of Joins Figure 4 Dissemination progress

Figure 1

 1 Figure 1 Effect of Tuples Arrival Rate Figure 2 Effect of Query Arrival Rate

Figure 5

 5 Figure 5 Speed of dissemination Figure 6 Message Overhead


Thomas Bréhard

Emmanuel Duflos
email: emmanuel.duflos@ec-lille.fr

Philippe Vanheeghe

Pierre-Arnaud Coquelin

Optimal Policies Search for Sensor Management : Application to the ESA Radar

Keywords: Sensor(s) Management, Partially Observable Markov Decision Process, Stochastic Gradient Estimation, AESA Radar. Sensor(s) Management Special Session

This paper introduces a new approach to solve sensor management problems. Classically sensor management problems can be well formalized as Partially-Observed Markov Decision Processes (POMPD). The original approach developped here consists in deriving the optimal parameterized policy based on stochastic gradient estimation. We assume in this work that it is possible to learn the optimal policy off-line (in simulation ) using models of the environement and of the sensor(s). The learned policy can then be used to manage the sensor(s). In order to approximate the gradient in a stochastic context, we introduce a new method to approximate the gradient, based on Infinitesimal Approximation (IPA). The effectiveness of this general framework is illustrated by the managing of an Electronically Scanned Array Radar.

I. INTRODUCTION

Years after years the complexity and the performances of many sensors have increased leading to more and more complex sensor(s)-based systems which supply the decision centers with an increasing amount of data. The number, the types and the agility of sensors along with the increased quality of data far outstrip the ability of a human to manage them: it is often difficult to compare how much information can be gained by way of a given management scheme [START_REF] Kalandros | Tutorial on multisensor management and fusion algorithms for target tracking[END_REF]. It results from this the necessity to derive unmanned sensing platforms that have the capacity to adapt to their environment [START_REF] Ji | Nonmyopic multiaspect sensing with partially observable markov decision processes[END_REF]. This problem is often refered as the Sensor(s) Management Problem. In more simple situations, the operational context may lead to works on sensor(s) management like in the radar -infrared sensor case [START_REF] Blackman | Design and Analysis of Modern Tracking Systems[END_REF]. A general definition of this problem could then be : sensor management is the effective use of available sensing and database capabilities to meet the mission goals. Many applications deal with military applications,a classical one being to detect, to track and tp identify smart targets (a smart target can change its way of moving or its way of sensing when it detects it is under analysis) with several sensors. The questions are then the following at each time: how must we group the sensors, how long, in which direction, and with which functioning mode? The increasing complexity of the targets to be detected, tracked and identified, makes the management even more difficult and led to the development of researches on the definition of an optimal sensor management scheme in which the targets and the sensors are treated altogether in a complex dynamic system [START_REF] Doucet | Particle filtering for multitarget tracking and sensor management[END_REF].

Sensor Management has become very popular this last years and many approaches can be found in the litterature. In [START_REF] Duflos | Time allocation of a set of radars in a multitarget environment[END_REF] and [START_REF] Huguerre | An optimal detection strategy for esa radars[END_REF] the authors use a the modelling of the detection process of an Electronically Scanned Array (ESA) Radar to propose management scheme during the detection step. In [START_REF] Kastella | Discrimination gain to optimize detection and classification[END_REF]- [START_REF] Khodayari-Rostamabad | Information theoric enumeration and tracking of multiple sources[END_REF] an information-based approach is use to manage a set of sensors. From a theorical point of view the sensor management can be modelled as a Partially Observable Markov Decision Process (POMDP) [START_REF] Kreucher | Adaptive multimodality sensor scheduling for detection and tracking of smart targets[END_REF]- [START_REF] Krishnamurthy | Structured threshold policies for dynamic sensor scheduling -a partially observed markov decision process approach[END_REF]. Whatever the underlying application, the sensor management problem consists in choosing at each time t an action A t within the set A of available actions. The choice of A t is generally based on the density state vector X t describing the environment of the system and variables of the system itself. It is generally assumed that the state or at least a part of this state is Markovian. Moreover in most of the applications, we only have access to a partial information of the state and X t must be estimated from the measurements {Y s } 1≤s≤t . This estimation process is often derived within a Bayesian framework where we use statedynamics and observation models such as:

X t+1 = F (X t , A t , N t ) (1)
Y t = H(X t , W t ) (2)
where N t , W t , F and H respectively stands for the state noise, the measurements noise, the state-dynamics and the measurement function. F and H are generally time varying functions. The control problem consists in finding the scheduling policy π i.e. select A t given the past and the possible futures. However, this control problem may have a theorical solution, it is generally untractable in practice. However few works propose optimal solution in the frame of POMDPs like [START_REF] Krishnamurthy | Structured threshold policies for dynamic sensor scheduling -a partially observed markov decision process approach[END_REF]. Beside, several works have been carried out to find sub-optimal policies like for instance myopic policies. Reinforcement Learning and Q-Learning have also been used to propose a solution ( [START_REF] Sutton | Time-derivative models of pavlovian reinforcement[END_REF], [START_REF] Kreucher | Non-myopic approaches to scheduling agile sensors for multitarget detection, tracking, and identification[END_REF]).

We propose in this paper to look for a policy within a class of parametrized policy π θ and to learn it which means learn the optimal value of θ. Funding our work on the approach described in [START_REF] Coquelin | Particle filter -based policy gradient in pomdps[END_REF] we assume that it is possible to learn this policy in simulation using models of the overall system. Once the optimal parameter has been found it is used to manage the sensor(s). The frame of this work being the detection and localization of targets, we show in the last part of this paper how it may be applied the the management of an ESA radar.

The section II described the modelling of a sensor management problem using a POMDP approach. In the section III we derive the algorithm to learn the parameter of the policy.In section IV we show how this method may be used for the tasking of an ESA radar. Finally section V exhibits firts simulations results.

II. MODELLING

A. POMDP Modelling

Let us consider three measurable continuous spaces denoted by X , A and Y, X is called the state space, Y the observation space and A the action space. We call M(X ) the set of all the measures defined on X . A Partially-Observable Decision Proceess is defined by a state process (X t ) t≥0 ∈ X , an observation process (Y t ) t≥1 ∈ Y and a set of actions (A t ) t≥1 ∈ A. In these definitions t stands usually for the time. The state process is an homogeneous Markov chain with initial probability measure µ(dx 0 ) ∈ M(X ) and with Markov transition kernel K(dx t+1 |x t ) ( [START_REF] Coquelin | Numerical methods for sensitivity analysis of feynman-kac models[END_REF]):

∀t ≥ 0, X t+1 ∼ K(•|X t ) (3)
X 0 ∼ µ (4)
(Y t ) t≥1 is called the observation is linked with the state process by the conditional probability measure:

P(Y t ∈ dy t |X t = x t ) = g(x t , y t ) dy t (5)
where g : X × Y → [0, 1] is the marginal density function of Y t given X t . In a general way, the state process evolves continuously with respect to time t whereas the observations are made at sampled time t n . A new observation is used to derive a new action. We will therefore consider in the following the processes (X t ) t≥0 , (Y n ) n∈N , (A n ) n∈N where n stands for the index of the observation. We also assume that there exists two generative functions F µ : U → X and F : X × U → X , where (U, σ(U ), ν) is a probability space, such that for any measurable test function f defined over X we have:

X f (x t )K(dx t |x t-1 ) = f (F (x t-1 , u))ν(du) (6)
and

X f (x 0 )µ(dx 0 ) = f (F µ (u))ν(du). (7)
In many practical situations, U = [0, 1] nU , and u is a n Uuple of pseudo random numbers generated by a computer. For sake of simplicity, we adopt the notations K(dx 0 |x -1 ) µ(dx 0 ) and F (x -1 , u) F µ (u). Under this framework, the Markov Chain (X t ) t≥0 is fully specified by the following dynamical equation:

X t+1 = F (X t , U t ), U t i.i.d. ∼ ν (8)
The observation process (Y n ) n∈N , defined on the measurable space (Y, σ(Y )), is actually linked with the state process by the following conditional probability measure

P(Y n ∈ dy n |X tn = x tn , A n ) = g(y n , x tn , A n )λ(dy t ) (9)
where A n ∈ A is defined on the measurable space (A, σ(A)) and λ is a fixed probability measure on (Y, σ(Y )). As we assume that observations are conditionally independent given the state process we cab write ∀1 ≤ i, j ≤ t, i = j:

P(Y i ∈ dy i , Y j ∈ dy j |X 0:t , A i , A j ) = P(Y i ∈ dy i |X 0:t , A i )P(Y j ∈ dy j |X 0:t , A j ) (10)
where we have adopted the usual notation z i:j = (z k ) i≤k≤j .

B. Filtering distribution in a Partially-Observable Markov Decision Process

Given a sequence of action A 1:n and a sample trajectory of the observation process y 1:n and indices

{n 1 , n 2 , t 1 , t 2 } such that 1 ≤ n 1 ≤ n 2 ≤ n and 0 ≤ t 1 ≤ t n1 ≤ t n2 ≤ t 2 ≤ t n ,
we define, using the , the posterior probability distribution M t1:t2|n1:n2 (dx t1:t2 ) by ( [START_REF]Feynman-Kac Formulae Genealogical and Interacting Particle Systems with Applications[END_REF]):

P(X t1:t2 ∈ dx t1:t2 |Y n1:n2 = y n1:n2 , A n1:n2 ) (11)
Using the Feynman-Kac framework, the probabilty 11 can be written:

t2 t=t1 K(dx t |x t-1 ) n2 j=n1 G tj (x tj ) X t 2 -t 1 t2 t=t1 K(dx t |x t-1 ) n2 j=n1 G tj (x tj ) , (12)
where for simplicity's sake, G tn (x tn ) g(y n , x tn , A n ) and G 0 (x 0 ) 0. One of the main interest here is to estimate the state at time t from noisy observations y 1:nt with n t the index of the last observation just before time t. From a bayesian point of view this information is completely contained in the so-called filtering distribution M t:t|1:nt . In the following, the filtering distribution will simply be denoted as M t .

C. Numerical methods for estimating the filtering distribution

Given a measurable test function f : X → R, we want to evaluate

M t (f ) = E[f (X t )|Y 1:nt = y 1:nt , A 1:nt ] (13)
which is equal, using the Feynman Kac framework, to:

E[f (X t ) nt j=1 G tj (X tj )] E[ nt j=1 G tj (X tj )] (14)
In general, it is impossible to find M t (f ) exactly except for simple cases such as linear/gaussian (using Kalman filter) or for finite state space Hidden Markov Models. In the general dynamics, continuous space case considered here, possible numerical methods for computing M t (f ) include the Extended Kalman filter, quantization methods, Markov Chain Monte Carlo methods and Sequential Monte Carlo methods (SMC), also called particle filtering. The basic SMC method, called Bootstrap Filter, approximates M t (f ) by an empirical distribution [START_REF] Doucet | On Sequential Monte Carlo Sampling Methods for Bayesian Filtering[END_REF]). It consists in a sequence of transition/selection steps: at time t, given observation y t ( [START_REF] Coquelin | Particle filter -based policy gradient in pomdps[END_REF]):

M N t (f ) = 1 N N i=1 f (x N i ) made of N so-called particles ( [
• Transition step: (also called importance sampling or mutation) a successor particles population x 1:N t is generated according to the state dynamics from the previous population x 1:N t-1 . The (importance sampling) weights w Resampling is used to avoid the problem of degeneracy of the algorithm, i.e. that most of the weights decreases to zero. It consists in selecting new particle positions such as to preserve a consistency property :

N i=1 w i t φ( x i t ) = E[ 1 N N i=1 φ(x i t )] (15)
The simplest version introduced in [START_REF] Gordon | Novel approach to nonlinear and non-gaussian bayesian state[END_REF] consists in choosing the selection indices k 1:N t by an independent sampling from the set 1 : N according to a multinomial distribution with parameters w 1:N t , i.e. P(k i t = j) = w j t , for all 1 ≤ i ≤ N . The idea is to replicate the particles in proportion to their weights. The reader can find some convergence results of

M N t (f ) to M t (f ) (e.
g. Law of Large Numbers or Central Limit Theorems) in [START_REF]Feynman-Kac Formulae Genealogical and Interacting Particle Systems with Applications[END_REF], but for our purpose we note that under weak conditions on the test function and on the HMM dynamics, we have the asymptotic consistency property in probability, i.e. lim N →∞ M N t (f )

P = M t (f ).

III. POLICY LEARNING ALGORITHM A. Optimal Parameterized Policy for Partially-Observable Markov Decision Process

Let R t be a real value reward function

R t R(X t , M t (f )) . ( 16
)
The goal is to find a policy

π : A n × Y n → A (17)
that maximizes the criterion performance :

J π = T 0 E[R t ]dt ( 18
)
where T is the duration of the scenario. Designing in practice policies that depend on the whole trajectory of the past observations/actions is unrealistic. It has been proved that the class of stationary policies that depend on the filtering distribution conditionally to past observations/actions M t contains the optimal policy. In general the filtering distribution is an infinite dimensional object, and it cannot be represented in a computer and so is the policy. We therefore propose to look for the optimal policy in a class of parameterized policies (π α ) α∈Γ that depend on a statistic of the filtering distribution :

A n+1 = π α (M tn (f )) ( 19
)
where f is any test function. As the policy π is parameterized by α, the performance criterion now depends only on α.

Thus we can maximize it by achieving a stochastic gradient ascent with respect to α :

α k+1 = α k + η k ∇J α k , k ≥ 0 (20)
where ∇J α k denotes the gradient of J α k w.r.t α k . By convention ∇J α k is column vector whose i-th component is the partial derivative with respect to α i . (η k ) k≥0 is a nonincreasing positive sequence tending to zero. We present in the two following subsection a possible approach to estimate ∇J α k based on Infinitesimal Perturbation Analysis (IPA).

B. Infinitesimal Perturbation Analysis for gradient estimation

We assume that we can write the following equality at each k:

∇J α = T 0 ∇ α E[R t ]dt (21)
Proposition 1: We have the following decomposition of the gradient

∇ α E[R t ] = E[M t (f S t )∇ Mt(f ) R t ] -E[M t (f )M t (S t )∇ Mt(f ) R t ] + E[R t S t ] (22)
where

S t = pt j=1 ∇ α G tj (X tj ) G tj (X tj ) (23)
Proof: First let us rewrite ∇ α E[R t ] as following:

∇ α E[R t ] = ∇ α X t ×Y n t R t U t V t nt j=1 λ(dy j ) (24)
where:

U t = t i=0 K(dx i |x i-1 ) , V t = nt j=1 G tj (x tj ) . ( 25
)
Remarking that only R t and V t depends on α so that we obtain

∇ α V t = S t V t , ∇ α R t = ∇ α M t (f )∇ Mt(f ) R t (26)
where S t is given by eq.( 23). Incorporating [START_REF] Vilmorin | Contribution à la grstion optimale de capteurs: application à la tenue de situations aériennes[END_REF] in (26), we obtain

∇ α E[R t ] = E[∇ α M t (f )∇ Mt(f ) R t ] . + E[R t S t ] . ( 27
)
Now using one more time (26), we have

∇ α M t (f ) = ∇ α E f (X t ) V t E[V t ] = E f (X t ) ∇ α V t E[V t ] -E f (X t ) V t E[∇ α V t ] E[V t ] 2 = E f (X t )S t V t E[V t ] -M t S t E V t E[V t ] = M t (f S t ) -M t (f )M t (S t ) (28)
so that we obtain [START_REF] Li | A Survey of Maneuvering Target Tracking Part I: Dynamics Models[END_REF] by incorporating (28) in ( 27).

We can deduce directly Algorithm 1 from [START_REF] Li | A Survey of Maneuvering Target Tracking Part I: Dynamics Models[END_REF]. It is important to note that we must deal with two time-scales. This first and the shorter one allows to simulate the continuous state X t . On the contrary the observation and action process are updated only each time we get a new observation. These specific time is denoted t n in the algorithm. That is the eason while there is an alternative to update the variables S t and w(i) t-1 . A new action A n is also calculated each t n as already explained above. One can also be surprised to calculate R(X t , M t (f )) using the sampled value of X t . To well understand this algorithm we must remind that the learning is made off-line using a simulated process. It is therefore possible to use the real value of X t in this case.

IV. APPLICATION TO THE ESA RADAR

The ESA is an agile beam radar which means that it is able to point its beam in any direction of the environnement almost instantaneously without inertia. However, the targets in the environement are detected w.r.t a probability of detection which depends on the direction of the beam and the time of observation in this direction. In the following, we precise first the nature of an action, then the influence of the action onto the probability of detection and finally the nature of the observations.

Definition of the action: The main property of an ESA is that it can point its beam without mechanically adjusting the antenna. An ESA radar provides measurements in a direction θ. We note δ, the time of observation in this direction. In this work the the n-th action is :

A n = θ n δ n T (29) Algorithm 1 Policy Gradient in POMDP via IPA Initialize α 0 ∈ Γ for k = 1 to ∞ do for t = 1 to T do Sample u t ∼ ν Set x t = F (x t-1 , u t ), If t = t n , sample y n ∼ g(., x t , a n )λ(.) Set s t = s t-1 + ∂g ∂α (xt,yn,an) g(xt,yn,an) if t = t n s t-1 else Set ∀i ∈ {1 . . . , I} x(i) t = F (x (i) t-1 , a t-1 , u (i) t ) where u (i) iid ∼ ν s(i) t = s i t-1 + ∂g ∂α (x i t ,yn,an) g(x (i) t ,yn,an) if t = t n s i t-1 else w(i) t =    g(x (i) t ,yn,an) w(i) t-1 P j g(x (j) t ,yn,an) w(j) t-1 if t = t n w(i) t-1 else Set (x (i) t , s (i) t ) i∈{1,...,I} = (x (i) t , s(i) t ) i∈{k1,...,kI } , k 1:I are selection indices associated to ( w(i) ) i∈{1,...,I} , m t (f ) = 1 I i f (x (i) t ), m t (s t ) = 1 I i s (i) , m t (f s t ) = 1 I i f (x (i) t )s (i) t , a n+1 = π α k (m t ) if t = t n r t = R(x t , m t (f )) ∇r t = (m t (f s t ) -m t (f )m t (s t )) ∂R ∂mt(f ) (x t , m t (f )) + r t s t ∇J α k = ∇J α k + ∇r t end for α k+1 = α k + η k ∇J α k end for with θ n ∈ [-π 2 , π 2 ] , δ n ∈ R + ∀n ≥ 0 . ( 30
)
This is a simple possible action. One could increase the number of componenets of an action by adding the emitted frequency for instance. The action does not influence directly the observation produced by the ESA but the probability of detection of a target.

The probability of detection P d : It refers to the probability to detect a target and therefore to the probability to obtain an estimation of the state of a target p at time t n denoted X tn,p with action A n . In this work, X tn,p is composed of the localisation and velocity components of the target p at time t n in the x-y plane: X tn,p = rx tn,p ry tn,p vx tn,p ry tn,p T (31) where the subscript T stands for matrix transpose. The terms rx tn,p and ry tn,p refers here to the position and vx tn,p and vy tn,p the velocity of target p at time t n . We also denote D n,p the random variable which takes values 1 if the radar produces a detection (and therefore an estimation) for target p and 0 else :

D n = D n,1 . . . D n,P T . (32)
As said previously, this probability also depends on the time of observation δ n . Aerial targets being considering here, the reflectivity of a target can be modelled using a Swerling I model [START_REF] Curry | Radar System Performance Modeling[END_REF]. We then have the following relation between the probability of detection and the probability of false alarm P f a (i.e. the probability that the radar produce a detection knowing that there is no target) ( [START_REF] Duflos | Time allocation of a set of radars in a multitarget environment[END_REF], [START_REF] Wintenby | Resource allocation in airborn surveillance radar[END_REF]):

P d (x tn,p , A n ) = P 1 1+ρ(x tn ,p ,An ) f a (33)
where ρ(x tn,p , A n ) is the target signal-to-noise ratio. In the case of an ESA radar, it is equal to :

ρ(x tn,p , A n ) = αδ n cos 2 θ n r 4 tn,p e -(β tn ,p -θn ) 2 2B 2 (34)
where r tn,p is the target range and β tn,p the azimuth associated to target p at instant time t n . α is a coefficient which includes all the parameters of the sensor and B is the beamwidth of the radar. It is reminded in Appendix A how the equations 33 and 34 may be derived. If we make the assumption that all the detections are independant, we can write :

P(D n = d n |X tn = x tn , A n ) = P p P(D n,p = d n,p |X tn,p = x tn,p , A n ) (35)
where

P(D n,p = d n,p |X tn,p = x tn,p , A n ) = P d (x tn,p , A n )δ dn,p=1 + (1 -P d (x tn,p , A n ))δ dn,p=0 (36)
Observation equation: At instant time t n , the radar produces a raw observation Y n composed of P measurements :

Y n = Y n,1 . . . Y n,P T . (37)
where Y n,p is the observation related to target of state value x tn,p obtained with action A n (we do not consider here the problem of measurement-target association). Moveover, we assume that the number of targets P is known. Each of these measurements has the following formulation :

Y n,p = r n,p β n,p ṙn,p T (38)
where r n,p , β n,p , ṙn,p are range, azimuth and range rate. The equation observation can be written

P(Y n ∈ dy n |X tn = x tn , A n ) = (39)
The relation between the state and the raw observations is given by :

Y n,p = h tn (X tn,p ) + W n,p (45)
with h tn (x tn,p ) equals to: and W n,p a gaussian noise the covariance matrix of which is given by :

      (rx
Σ y = diag(σ 2 r , σ 2 β , σ 2 ṙ ) . ( 47
)
State equation: First let us introduce the definition of the unknown state X t at time t and its evolution through time. X t,p is the state of the target p. It has been defined above. Let P be the known number of targets in the space under analysis at time t. X t has the following form: .

X t = X t,1 . . . X t,P T (48)
Based on [START_REF] Li | A Survey of Maneuvering Target Tracking Part I: Dynamics Models[END_REF] works, we classically assume that all the targets follow a nearly constant velocity model. We use a discretized version of this model ( [START_REF] Cadre | Bearings-only tracking for maneuvering sources[END_REF]) :

X t,p = F (X t-1,p , U t ) where U t ∼ N 0, σ 2 Q (49)
where

F =     1 0 β 0 0 1 0 β 0 0 1 0 0 0 0 1     and Q =      β 3 3 0 β 2 2 0 0 β 3 3 0 β 2 2 β 2 2 0 β 0 0 β 2 2 0 β      . (50)
V. SIMULATIONS

VI. CONCLUSION APPENDIX A

We show in this Appendix how the probability of detection is derived. First, the radar transmits a pulse expressed as follows

s(t) = α(t) cos(w c t) (51) = Re{α(t)e jwct } ( 52
)
where α(t) is the envelope also called the transmitted pulse and w c the carrier frequency. This pulse is modified by the process of reflection. A target is modelled as a set of elementary reflectors, each reflecting: time delayed, Doppler shift, Phase shift and attenuated version of the transmitted signal. We usually assume that the reflection process is linear and frequency independent within the bandwidth of the transmitted pulse. The return signal has the following formulation:

s r (t) = G i α(t -τ i )g i e j(wc(t-τi+ 2 ṙi c t)+θi) + n(t) (53)
where

• g i is the radar cross section associated to reflector i,

• θ i is the phase shift associated to reflector i,

• ṙi is the radial velocity between the antenna and the object (Doppler frequency shift), • G: others losses heavily range dependent due to spatial spreading of energy, • n(t) is a thermal noise of the receiver such that Re{n(t)}, Im{n(t)} ∼ N (0, σ 2 n ). We make the following approximations:

ṙi ≈ ṙ α(t -τ i ) ≈ α(t -τ ) ( 54
)
where ṙ is the mean radial velocity of the target τ is the mean time delay of the target. Using these approximations, the return signal can be rewritten as follows:

s r (t) = α(t -τ )Ge jwD t b + n(t) (55)
where

w D = w c (1 + 2 ṙi c ) b = i g i e j(-wcτi+θi) . (56)
The fluctuations of b are known and modelled using Swerling 1 model [START_REF] Curry | Radar System Performance Modeling[END_REF]. There are differents models availables (Swerling 1, 2, 3,...) corresponding to different types of targets. Swerling 1 given below is convenient for aircrafts. We can then write :

Re{b}, Im{b} ∼ N (0, σ 2 RCS ) . (57)
This modelling of b assumes that the phase shifts θ i are independent and uniformly distributed and the magnitudes g i are identically distributed. If the number of reflector is large, the central limit theorem gives that b is a complexvalued Gaussian random variable centered at the origin. Now, a matching filter is applied to our return signal

s m (t) = +∞ -∞ s r (t)h(s)ds (58)
where h(t) is a shifted, scaled and reversed copy of s r (t)

h(s) = α(δ -t)e -jwD (δ-t) . (59)
We choose t = δ +τ which yields the best signal to noise ratio where δ is the length of the transmitted pulse. The probability of detection is based on quantity |s m (δ + τ )| 2 . We can show that

s m (δ + τ ) = Ge jwD τ b + +∞ -∞ n(δ + τ -s)h(s)ds . ( 60
)
One can remark that s m (δ + τ ) is the sum of two complexvalue Gaussian variables. We look at the following statistic is called the Signal-to-Noise Ration noted ρ. This SNR is related to the parameters of the system and the target. The classical radar equation is given by the following formula ( [START_REF] Wintenby | Resource allocation in airborn surveillance radar[END_REF]):

ρ = P t G t G r λ 2 σ (4π) 3 r 4 ( 67
)
where P t is the energy of the transmitted pulse, G t is the gain of the transmitted antenna, G r is the gain of the received antenna, σ is the radar cross section (for an aircraft between 0.1 and 1 m 2 ), r is the target range, γ is the system noise temperature and L is a general loss term. However, the above formula does not take into account for the sake of simplicity the losses due to atmospheric attenuation and to the imperfection of the radar. Thus , extra terms must be added :

ρ = P t G t G r λ 2 σ (4π) 3 kbLγr 4 ( 68
)
where b is the receiver noise bandwith (generally consider equal to the signal bandwidth so that b = 1 δt ), k is Boltzmann's constant, γ is the temperature of the system and L some losses. Moreover, the gain reduces with the deviation of the beam from the antenna normal in an array antenna.

G t = G 0 cos α (θ t ) , ( 69
) G r = G 0 cos α (θ t ) (70)
where G 0 is the gain of the antenna. In [START_REF] Vilmorin | Contribution à la grstion optimale de capteurs: application à la tenue de situations aériennes[END_REF], α = 2, in [START_REF] Wintenby | Resource allocation in airborn surveillance radar[END_REF], α = 2.7. According [START_REF] Keuk | On Phased-Array Radar Tracking and Parameter Control[END_REF], there is also a beam loss because the radar beam is not pointing directly so that the radar equation is:

ρ = P t G 2 0 λ 2 σδ t cos 2 (θ t ) (4π) 3 kLγr 4 e -(θ t -β t ) 2 2B 2 (71)
where is B is the beamwidth.

  g(e x j t ,yt) are evaluated. • Selection step: Resample (with replacement) N particles x 1:N t from the set x 1:N t according to the weights w 1:N

H 1 :H 0 P

 10 data generated by signal + noise H 0 : data generated by noise (:p Λ (x) = e -x(64)Then, we derive the probability of detection and false alarm.   P d = +∞ γ p Λ (x|H 1 is true) = e f a = +∞ γ p Λ (x|H 0 is true ) = e -γ

  ∈ dy n,p |X tn,p = x tn,p , A n ) (41) = g(y n,p , x tn,p , A n )λ(dy n,p )(42) g(y n,p , x tn,p , A n ) = N (h t (x tn,p ), Σ y ) P d (x tn,p , A n ) 1 -P d (x tn,p , A n ) ) = λ cont (dy n,p ) + λ disc (dy n,p )

	P p P(Y n,p T P(Y n,p ∈ dy n,p |X tn,p = x tn,p , A n ) (40) where and λ(dy n,p	(43)

  tn,prx obs tn ) 2 + (ry tn,pry obs tn ) 2

			
	(rxt n ,p-rx	atan tn )(vxt n ,p -vx obs ryt n ,p -ry obs tn rxt n ,p -rx obs tn obs tn )+(ryt n,p -ry obs tn )(vyt n ,p -vy obs tn ) (rxt n,p -rx obs tn ) 2 +(ryt n,p -ry obs tn ) 2 √	    	(46)


Bartłomiej Nowak
email: bartlomiej.nowak@pwr.edu.pl

Michel Grabisch
email: michel.grabisch@univ-paris1.fr

Katarzyna Sznajd-Weron
email: katarzyna.weron@pwr.edu.pl

The threshold model with anticonformity under random sequential updating

We study an asymmetric version of the threshold model of binary decision making with anticonformity under asynchronous update mode that mimics continuous time. We analyze this model on a complete graph using three different approaches: the mean-field approximation, Monte Carlo simulation, and the Markov chain approach. The latter approach yields analytical results for arbitrarily small systems, in contrast to the mean-field approach, which is strictly correct only for an infinite system. We show that for sufficiently large systems, all three approaches produce the same results, as expected. We consider two cases: (1) homogeneous, in which all agents have the same tolerance threshold, and (2) heterogeneous, in which thresholds are given by a beta distribution parameterized by two positive shape parameters α and β. The heterogeneous case can be treated as a generalized model that reduces to a homogeneous model in special cases. We show that particularly interesting behaviors, including social hysteresis and critical mass reported in innovation diffusion, arise only for values of α and β that yield the shape of the distribution observed in reality.

I. INTRODUCTION

Within the broad class of two-state dynamics [START_REF] Gleeson | Binary-state dynamics on complex networks: Pair approximation and beyond[END_REF], threshold models are particularly useful for describing various social and economic phenomena [START_REF] Granovetter | Threshold models of collective behavior[END_REF][START_REF] Watts | A simple model of global cascades on random networks[END_REF][START_REF] Breer | Models of tolerant threshold behavior (from T. Schelling to M. Granovetter[END_REF]. As other binary-state opinion dynamics [START_REF] Jędrzejewski | Statistical Physics Of Opinion Formation: Is it a SPOOF?[END_REF], the threshold model describes the social influence in decision making for the choice between precisely two alternatives, often denoted by 1 (agree, adopt, be active, etc.) and 0 (disagree, refuse, be inactive, etc.). Although a binary decision framework seems to be oversimplified, it is relevant for surprisingly many complex problems [START_REF] Watts | A simple model of global cascades on random networks[END_REF].

In the original threshold models of collective behavior proposed by Schelling [START_REF] Schelling | Micromotives and macrobehavior[END_REF] and Granovetter [START_REF] Granovetter | Threshold models of collective behavior[END_REF], an agent takes action 1 if the proportion of his neighbors in state 1 exceeds some threshold, otherwise action 0 is taken. It means that an agent in state 1 may return to state 0, because not enough neighbors are active. On the other hand, in many other threshold models, the transition from state 1 to state 0 is forbidden [START_REF] Watts | A simple model of global cascades on random networks[END_REF][START_REF] Dodds | Universal behavior in a generalized model of contagion[END_REF][START_REF] Juul | Hipsters on networks: How a minority group of individuals can lead to an antiestablishment majority[END_REF].

Here, we will use the original formulation, in which a transition from 1 to 0 is possible, as in [START_REF] Granovetter | Threshold models of collective behavior[END_REF][START_REF] Lee | Social contagion with degreedependent thresholds[END_REF][START_REF] Grabisch | Anti-conformism in the threshold model of collective behavior[END_REF], but additionally in the presence of anticonformity. Such a model has been already studied from a mathematical point of view under the synchronous update mode [START_REF] Grabisch | Anti-conformism in the threshold model of collective behavior[END_REF], considering a complete network or a random neighborhood. The study focused on finding absorbing classes, cycles, etc. In this paper, we investigate the same model but under random sequential updating (asynchronous updating), which mimics continuous time, and restricting to the complete network. Quick analysis shows that the two models behave very differently. Under the complete network assumption, the synchronous updating is a deterministic model and causes the appearance of cycling phenomena, in a way which is very dependent on the parameters chosen. On the contrary, asynchronous updating is probabilistic and prevents the appearance of cycles. We show that it behaves similarly to a random-walk process. The appropriate tools in this situation are the phase transitions and phase diagrams under the mean-field approximation, which are typical for statistical physics of opinion formation [START_REF] Castellano | Statistical physics of social dynamics[END_REF][START_REF] Galam | Two-dimensional Ising transition through a technique from two-state opiniondynamics models[END_REF][START_REF] Raducha | Statistical mechanics of a coevolving spin system[END_REF][START_REF] Calvelli | Phase transitions and universality in the Sznajd model with anticonformity[END_REF][START_REF] Vieira | Pair approximation for the noisy threshold q-voter model[END_REF], and also the Markov chain approach. It should be noted that phase transition diagrams cannot be used in the case of cycling, as only steady-states can be studied.

We study the model on a complete graph since in this case the mean-field method allows us to obtain rigorous results. Independently, we conducted Monte Carlo simulations to validate the theoretical approach. Finally, we present a Markov chain approach, which allows us not only to obtain results for arbitrary small systems but also to derive the stationary distribution of visited states.

II. MODEL

We consider a society of n agents placed at the vertices of an arbitrary graph G = (N, E), where N = {1, ..., n} is a set of vertices (agents) and E is the set of undirected edges. Each agent i has a set of neighbors K i = {j ∈ N : {i, j} ∈ E}, and the cardinality of this set |K i | = k i is the degree of agent i. Here we assume that each agent belongs to its own neighborhood to avoid the unrealistic effect that the agent does not consider its own state at all. However, this assumption is only relevant when the neighborhood is small, which is obviously not the case for an infinite complete graph. As in many other models, an agent can be in one of two alternative states: 1 (agree, adopt, be active, etc.) or 0 (disagree, refuse, be inactive, etc.). Following [START_REF] Grabisch | Anti-conformism in the threshold model of collective behavior[END_REF], we use the term 'active' for agents in state 1, and 'inactive' for agents in state 0, and denote by a i (t) the state (action) taken by agent i at time t.

We consider two types of social response, anticonformity and conformity, which occur with complementary probabilities p and 1p respectively. In both cases, an agent can change its state if the ratio of active neighbors is above its tolerance threshold r i ∈ [0, 1]. The threshold r i of each agent is the realization of the random variable R with arbitrary distribution function F R (r) and does not change over time. In case of conformity, an agent follows the others, whereas in case of anticonformity he takes an opposite state to others. Therefore, the dynamics of the agent's state in the case of conformity can be written as [START_REF] Grabisch | Anti-conformism in the threshold model of collective behavior[END_REF]:

a i (t + ∆t) =    1, if 1 ki j∈Ki a j (t) ⩾ r i 0, otherwise, (1)
whereas in case of anticonformity [START_REF] Grabisch | Anti-conformism in the threshold model of collective behavior[END_REF]:

a i (t + ∆t) =    0, if 1 ki j∈Ki a j (t) ⩾ r i 1, otherwise. (2)
In this paper, we use the random sequential update mode, which means that an elementary update consists of:

(i) random drawing of agent i from all n agents (ii) with probability p agent i anticonforms to the neighborhood, that is, takes action a i (t + ∆t) according to Eq. (

(iii) with complementary probability 1p agent i conforms to the neighborhood, that is, takes action a i (t + ∆t) according to Eq. (

(iv) time is updated:

t := t + ∆t
As usually, ∆t = 1/n which means that the time unit consists of n elementary updates, which corresponds to one Monte Carlo step (MCS).

III. TRANSITION PROBABILITIES

Since we limit our study to the complete graph, we can fully describe the state of the system using a single random variable:

c = n 1 n , (3)
where n 1 is the number of agents in the state 1 and thus c is the ratio of active agents. Therefore, there are n + 1 possible states of the system: 0, 1 n , 2 n , . . . , 1. Because we use the sequential (asynchronous) update mode, at most one agent can change its state at a time, and thus we can introduce the following transition probabilities:

γ + (c) = P r c(t + ∆t) = c(t) + 1 n , γ -(c) = P r c(t + ∆t) = c(t) - 1 n . (4)
For our model, the explicit form of these probabilities can be written, according to the algorithm described in the previous section, as follows:

γ + (c) = (1 -p)(1 -c)P r(R ⩽ c) + p(1 -c)P r(R > c), γ -(c) = (1 -p)cP r(R > c) + pcP r(R ⩽ c), (5)
where P r(R ⩽ c) is the probability that the concentration c of the active agents is greater than or equal to the threshold R of the considered agent. This probability is simply the value of the cumulative distribution function F R (r) at r = c. Similarly, P r(R > c) is the probability that the concentration of active voters does not exceed the threshold of considered agents, and thus it is equal to 1 -F R (c). Therefore, we obtain

γ + (c) = (1 -p)(1 -c)F R (c) + p(1 -c)(1 -F R (c)), γ -(c) = (1 -p)c(1 -F R (c)) + pcF R (c). (6)
As can be seen from Eq.( 4), the concentration of active agents c is a random variable. However, we can easily write the evolution equation for the expected value of c. Moreover, for n → ∞ we can assume that c is localized to the expectation value. Therefore, we can write [START_REF] Jędrzejewski | Statistical Physics Of Opinion Formation: Is it a SPOOF?[END_REF]:

dc dt = γ + (c) -γ -(c) ≡ f (c), (7)
where f (c) can be interpreted as an effective force acting on the system. Such a force will later allow us to introduce a potential that helps visualize the dynamics of the system [START_REF] Strogatz | Nonlinear Dynamics and Chaos: With Applications to Physics[END_REF]. As usually, we focus on the steady states, i.e. those for which

dc dt = 0. (8)
In the next two sections, we will use condition [START_REF] Juul | Hipsters on networks: How a minority group of individuals can lead to an antiestablishment majority[END_REF] to calculate the stationary concentration of active agents for two cases: (1) homogeneous, in which all agents have the same tolerance threshold (2) heterogeneous, in which the distribution of thresholds F R (r) is given by the beta distribution. We will compare the analytical results with the results of Monte Carlo simulations for the system of size n = 10 4 , averaged over 10 independent runs collected after 10 4 Monte Carlo steps. For the Monte Carlo simulations, two types of initial conditions will be used to reproduce all stable solutions of Eq. ( 8): (1) all agents initially active, which will be denoted by c(0) = 1 and

(2) all agents initially inactive, which will be denoted by c(0) = 0. 16). Symbols represent Monte Carlo simulations from two initial conditions indicated in the legend. Arrows indicate the flow in the system. The length of each arrow correspond to the value of velocity dc dt given by Eq. [START_REF] Grabisch | Anti-conformism in the threshold model of collective behavior[END_REF].

IV. ONE THRESHOLD

In this case, the random variable R takes one value for all agents in the system, that is, all voters have the same threshold r:

F R (c) = 1 {r⩽c} , 1 -F R (c) = 1 {r>c} , (9)
where 1 {r⩽c} = 1 when r ⩽ c and 0 otherwise. Inserting (9) to Eq.( 6) and then into Eq. ( 7) we obtain:

dc dt = (1 -p) (1 -c)1 {r⩽c} -c1 {r>c} +p (1 -c)1 {r>c} -c1 {r⩽c} (10)
From [START_REF] Grabisch | Anti-conformism in the threshold model of collective behavior[END_REF], we obtain several trivial fixed points:

(p = 0, c = 0) ∀r ̸ = 0 (11) (p = 1, c = 0) r = 0 (12) (p = 0, c = 1) ∀r ∈ [0, 1] (13)
(p = 0.5, c = 0.5) ∀r ∈ [0, 1] (14)
(p = 1 -r, c = r) ∀r ∈ [0, 1] (15)
The remaining solutions can be obtained by solving Eq. ( 8), which leads to

p = 1 {r⩽c} -c 1 {r⩽c} -1 {r>c} , (16)
what is equivalent to the following cases,

∀c < r p = c, ∀c ⩾ r p = 1 -c. (17)
From the above analysis, we do not obtain the steady state for any value of r ∈ [0, 1] if p ⩾ r for r > 0.5 or p > 1r for r ⩽ 0.5. However, from the evolution of Eq. ( 10), represented by solid lines in Fig. 1, and from Monte Carlo simulations, represented by symbols in Fig. 1, we see that the system approaches the state c = r for p > r. This raises the question of what the evolution of the system actually looks like. To answer this question we can provide at least three alternative approaches. We decided to present all three because we believe that there is value in showing the different possibilities of system analysis, especially given the interdisciplinary nature of the work. The composition of the authors themselves is also interdisciplinary, and each of us found a different method more compelling/intuitive. The first, which is the most basic technique used to analyze dynamical systems, consists of interpreting a differential equation as a vector field [START_REF] Strogatz | Nonlinear Dynamics and Chaos: With Applications to Physics[END_REF]. Within this graphical way of thinking, we draw the arrows representing the flow dc/dt in the plane of model parameters, as shown in Fig. 1. The second method is based directly on the transition probabilities γ + (c), γ -(c). If we split the transition probabilities into cases,

∀c < r γ + (c) = p(1 -c) ∧ γ -(c) = (1 -p)c, ∀c ⩾ r γ + (c) = (1 -p)(1 -c) ∧ γ -(c) = pc. (18)
we easily observe that they do not cross at any point, when ∀r > 0.5 p ⩾ r or ∀r ⩽ 0.5 p > 1r, see the fourth column of Fig. 2. This implies no steady state. On the other hand, for p < r transition probabilities γ -(c) and γ + (c) cross each other, as shown in the first three columns of Fig. 2, i.e., the steady state γ -(c) = γ + (c) exists.

The third method is based on the idea of potential V (c) [START_REF] Strogatz | Nonlinear Dynamics and Chaos: With Applications to Physics[END_REF]:

V (c) = -f (c)dc. ( 19
)
From Eq. [START_REF] Yang | Dynamical system model predicts when social learners impair collective performance[END_REF] we see that the effective force f (c) defined in Eq. ( 7) takes the following form:

f (c) = γ + (c) -γ -(c) = p -c, for c < r 1 -p -c, for c ⩾ r (20)
and thus: FIG. 2. Analysis of the steady states and the stability of the system for two values of threshold r = 0.6 (two first rows) and r = 0.5 (two last rows) and four values of p = 0.2 (first column), p = 0.4 (second column), p = 0.45 (third column), p = 0.8 (fourth column). In first and third row, solid lines represents values of γ + and dotted lines stands for γ -obtained with Eqs. [START_REF] Yang | Dynamical system model predicts when social learners impair collective performance[END_REF]. Potentials V (c) (second and fourth row) are obtained with Eqs. [START_REF] Seneta | Non-negative Matrices and Markov Chains[END_REF]. In all subplots, filled circles denote continuity, while empty circles denote lack of continuity at this point.

V (c) = -(p -c)dc = c 2 2 -cp, for c < r -(1 -p -c)dc = c 2 2 -c(1 -p), for c ⩾ r (21) 0 γ + (c), γ -(c) 0 γ + (c), γ -(c) 0 γ + (c), γ -(c) 0 γ + (c), γ - (c
0 c V (c) 0 γ + (c), γ -(c) 0 γ + (c), γ -(c) 0 γ + (c), γ -(c) 0 γ + (c), γ - (c
Using such an approach, we draw a ball sliding down the walls of a potential well [START_REF] Strogatz | Nonlinear Dynamics and Chaos: With Applications to Physics[END_REF], as shown in Fig. 2. It should be noted that potentials given by Eq. ( 21) are determined up to a constant. Therefore, the size of at the border = not unambiguously Here, we have assumed that both constants are equal to zero, but this does not influence the dynamics of the system since the system cannot jump from one to the other potential minima, as presented in Fig. 2. This is because the point c = r belongs only to one area and thus the dynamics from this point is uniquely defined no matter what jump is at c = r.

The steady states are the local extrema of V (c). From Eq. ( 21) we see that the potential has a discontinuity at c = r, which implies no maximum (unstable steady state). Still, at most two minima (stable steady states) are possible. In general, the number of minima, denoted by M (r, p), can be described as follows:

∀r > 0.5 M (r, p) =      2, for p ⩽ 1 -r 1, for 1 -r < p < r 0, for p ⩾ r, (22)
∀r ⩽ 0.5 M (r, p) =      2, for p < r 1, for r ⩽ p ⩽ 1 -r 0, for p > 1 -r. (23)
In conclusion, despite the lack of steady state in the case M (r, p) = 0 we can observe that the flow of the system is towards the point c = r. It reaches an asymptotic minimum at this point because from both the left and right boundaries, the system flow is towards this minimum. This explains the behavior shown in Fig. 1, which was at first incomprehensible and inspired the above analysis.

V. BETA DISTRIBUTION

In the previous section, we studied the homogeneous system, in which all agents had the same value of the tolerance threshold r. However, we can also consider more general distributions of thresholds, allowing for heterogeneity. The most useful are distributions whose support values r ∈ [0, 1] and show a variety of shapes. This is the case of the beta distribution with two parameters α and β, previously considered, for the models of tolerance without anticonformity [START_REF] Breer | Models of tolerant threshold behavior (from T. Schelling to M. Granovetter[END_REF]. It has a well-defined cumulative distribution function:

F R (r) = I r (α, β) = B(r, α, β) B(α, β) , (24)
where I r (α, β) is the regularized incomplete beta function, which can be defined in terms of the incomplete beta function B(r, α, β) and the complete beta function B(α, β). Inserting F R (r) given by Eq. ( 24) to ( 6) we obtain the transition probabilities γ + (c), γ -(c). Then inserting them to Eq. ( 7) we get

dc dt = (1 -p) [(1 -c)I c (α, β) -c(1 -I c (α, β))] +p [(1 -c)(1 -I c (α, β)) -cI c (α, β)] . (25)
Again, we can point out the obvious steady states (p = 0, c = 0), (p = 0, c = 1) for arbitrary values of α and β.

For c = 1/2 formula (25) boils down to

dc dt c= 1 2 = I 1 2 (α, β) - 1 2 (1 -2p), (26)
which has two roots. The first one p = 1/2 gives the fixed point (p = 1 2 , c = 1 2 ). The other root I 1 2 (α, β) = 1 2 exists if the beta distribution is symmetric around the value 1 2 . This happens for α = β, which leads to the conclusion that the value c = 1 2 is a fixed point for all values of p if α = β. For all remaining solutions the relation:

p = I c (α, β) -c 2I c (α, β) -1 . ( 27
)
0 1 0 1 β α β > α β < α β > α β < α β = α > 1 β = α < 1 β = α → ∞ β → ∞, β >> α β → ∞, β > α α → ∞, β < α α → ∞, β << α FIG. 3.
Phase diagram for the heterogeneous model with thresholds described by the beta distribution parametrized by two shape parameters α and β. Each inset shows representative behavior of c(p) for a given area of the phase diagram. Solid lines in the insets correspond to stable stationary states, whereas dashed lines correspond to unstable stationary states.

The information about the stability of the steady state is given by the sign of the derivative

dF dc = c α-1 (1 -c) β-1 Γ(α + β) Γ(α)Γ(β) (1 -2p) -1. ( 28
)
The state is stable if the above derivative is negative and unstable otherwise. The overall behavior of the model is summarized in Fig. 3. In the insets of this figure the dependence between the stationary value of c and parameter p is shown. Two shaded areas in Fig. 3 correspond to the situation in which at least one of the parameters α, β is smaller than 1. In this case, for p > 0 there is always only one steady state and c is monotonically increasing (β < α), monotonically decreasing (β > α) or constant (β = α) function of p.

Recalling the shape of the probability density function (PDF) of the beta distribution, we can draw some conclusions. If the PDF of the tolerance threshold is a monotonically decreasing function of the threshold r, then the concentration of active agents decreases with the probability of anticonformity p, and vice versa. If the PDF has the highest values at r = 0 and r = 1, being a convex function of r, then for all values of p > 0 the stationary value of active agents is 0.5.

The most complex behavior is seen if both shape parameters α, β are greater than 1 but not infinitely large, which corresponds to a uni-modal PDF, with zero probabilities at both end of the interval range, i.e., at r = 0 and r = 1. This case correspond to moderate tolerance [START_REF] Breer | Models of tolerant threshold behavior (from T. Schelling to M. Granovetter[END_REF], and it is a typical shape of the distribution of actual trait manifestations in behavior, as reported by psychologists [START_REF] Fleeson | The implications of Big Five standing for the distribution of trait manifestation in behavior: Fifteen experience-sampling studies and a meta-analysis[END_REF]. In such a case, the phase transitions appear, as shown in Fig. 3. As long as β = α, which corresponds to the symmetric PDF, there is a continuous phase transition between the phase in which one type of agent (active or inactive) dominates, and the symmetrical phase without the domination. The critical point, at which this transition occurs, can be calculated by solving the equation

dF dc c= 1 2 = 0, (29)
which gives:

p * 1 = 1 2 - Γ 2 (α) 2 3-2α Γ(2α) . ( 30
)
For α ̸ = β, as long as shape parameters are finite and at least one of them is larger than 1, we obtain an interesting behavior, with the jump at some value of p = p * and hysteresis, as shown in Fig. 3. This can be especially useful to describe the innovation diffusion. For example, if β > α then for the small value of p < p * there is possibility of high adoption if the initial fraction of adopted is above the critical mass. However, if the initial fraction of adopted is too low, i.e. below the critical mass, the innovation cannot spread in the society. Similar behavior has been recently reported for the completely different mathematical model of the collective decision-making with social learners for unequal merit options [START_REF] Yang | Dynamical system model predicts when social learners impair collective performance[END_REF].

It is worth noticing that for α, β → ∞ we can recover the solution for the model with one threshold, as shown in Fig. 3. We are able to do that by recalling the formula giving the mode of the beta distribution with α, β > 1:

m = α -1 α + β -2 . ( 31
)
While α, β → ∞, the beta distribution is a 1-point degenerate distribution with probability 1 at the midpoint m and 0 elsewhere. Thus, to obtain the case with the mode at the point m = r, i.e., recover the distribution for one threshold, parameters α and β should follow the formula

β = (1 -r)α -1 + 2r r (32) for α, β → ∞.
All results obtained analytically for beta distribution can be also obtained by Monte Carlo simulations, as shown in Fig. 4.

VI. MARKOV CHAIN APPROACH

Previously, we were assuming that the size of the system is infinite, i.e., n → ∞. However, such an assumption is not very realistic for social systems. Actually, social scientists are often interested in small systems. Therefore, in this section, we make analysis of the convergence of c in the long run using Markov chains for arbitrary small systems. The advantage of the Markov chain approach in the context of agent-based modeling of opinion dynamics has been already reported in [START_REF] Banisch | Agent based models and opinion dynamics as markov chains[END_REF].

Transition probabilities given by Eq. ( 4) allows us to write the (n+1)×(n+1) transition matrix, whose general term (i, j) indicates the probability of transition from state i to state j. Due to the asynchronous update mode, P is a tridiagonal row-stochastic matrix:

P =              γ 0 (0) γ + (0) 0 0 • • • 0 γ -1 n γ 0 1 n γ + 1 n 0 • • • 0 0 γ -2 n γ 0 2 n γ + 2 n • • • 0 0 • • • . . . . . . . . . 0 0 • • • 0 γ -n-1 n γ 0 n-1 n γ + n-1 n 0 • • • 0 0 γ -(1) γ 0 (1)              (33) with γ 0 (c) = 1 -γ + (c) -γ -(c).
This process is a random walk process. Its transition graph is strongly connected and aperiodic, hence P is a primitive matrix, i.e., the only absorbing class is the set of all states. This means that in the long run, the system at time t can be in any of the n + 1 states, and there is no stabilization [START_REF] Kemeny | Finite Markov chains[END_REF][START_REF] Seneta | Non-negative Matrices and Markov Chains[END_REF].

From Markov chain theory, the limit vector

π = [π(0), • • • , π(c), • • • , π(1)
] T giving the probability π(c) to be in state c in the long run is obtained as the left normalized eigenvector of P associated to eigenvalue 1, i.e., π is the solution of the linear system in variable z

(P T -I)z = 0, ( 34
)
1 T z = 1 (35)
with 1 = (1, . . . , 1) T . From now on, to avoid heavy no-tation, we denote γ + (k/n) by γ + (k), and similarly for γ -(k/n), π(k/n), etc. We obtain

P T -I =                  -γ + (0) γ -(1) 0 0 0 • • • 0 γ + (0) -γ -(1) -γ + (1) γ -(2) 0 0 • • • 0 0 . . . . . . . . . 0 • • • 0 0 • • • γ + (k -1) -γ -(k) -γ + (k) γ -(k + 1) • • • 0 0 • • • 0 . . . . . . . . . 0 0 • • • 0 0 γ + (n -2) -γ -(n -1) -γ + (n -1) γ -(n) 0 • • • 0 0 0 γ + (n -1) -γ -(n)                  (36)
Solving the system yields

π(0) = γ -(1) γ + (0) π(1) π(1) = γ -(2) γ + (1) π (2)
. . . = . . .

π(k) = γ -(k + 1) γ + (k) π(k + 1) (37)
. . . = . . .

π(n -1) = γ -(n) γ + (n -1) π(n).
This yields:

π(k) = γ + (k -1) γ -(k) γ + (k -2) γ -(k -1) • • • γ + (0) γ -(1) π(0) (k = 1, . . . , n).
(38) In the case of one threshold we are able to derive the above formulas analytically. Using Eqs. [START_REF] Yang | Dynamical system model predicts when social learners impair collective performance[END_REF] we obtain:

π(k) = (1 -p)(k + 1) p(n -k) π(k + 1) (k < rn -1) π(k) = k + 1 n -k π(k + 1) (rn -1 ⩽ k < rn) π(k) = p(k + 1) (1 -p)(n -k) π(k + 1) (k ⩾ rn).
Let us find when π(k) is increasing or decreasing. Supposing k < rn -1, we have:

(1 -p)(k + 1) p(n -k) ⩽ 1 ⇔ (1 -p)(k + 1) ⩽ p(n -k) ⇔ k ⩽ p(n + 1) -1.
When k ≥ rn, we obtain:

p(k + 1) (1 -p)(n -k) ⩽ 1 ⇔ k ⩽ n -p(n + 1).
Therefore, (i) For states below r, the peak is attained at

ĉ1 = k1 n , with k1 = ⌈p(n + 1)⌉ -1.
Observe that when n is large, this yields ĉ1 ≈ p.

(ii) For states above r, the peak is attained at

ĉ2 = k2 n , with k2 = n -⌊p(n + 1)⌋.
When n is large, we obtain ĉ2 ≈ 1p.

Depending on the relative positions of p and r, there can be one or two peaks, as summarized as follows:

• If r ⩽ p, r ⩽ 1 -p: peak at ĉ2 , • if p ⩽ r ⩽ 1 -p: two peaks at ĉ1 , ĉ2 , • if 1 -p ⩽ r ⩽ p: peak at ⌈rn⌉ n , • p ⩽ r, 1 -p ⩽ r: peak at ĉ1 .
In the case there are two peaks, i.e., p ≤ r ≤ 1p, let us find the relative heights of the peaks. From (38), we find, assuming rn ̸ ∈ N, Hence, assuming p(n + 1) ̸ ∈ N,

π(⌊rn⌋) = π k1 p 1 -p ⌊rn⌋-⌈p(n+1)⌉+1 × (n -⌊rn⌋ + 1) • • • (n -⌈p(n + 1)⌉ + 1) ⌊rn⌋ • • • ⌈p(n + 1)⌉ π(⌊rn⌋ + 1) = π k2 p 1 -p n-⌊p(n+1)⌋-⌊rn⌋-1 × (⌊rn⌋ + 2) • • • (n -⌊p(n + 1)⌋) (n -⌊rn⌋ -1) • • • (⌊p(n + 1)⌋ + 1) π(⌊rn⌋ + 1) = π(⌊rn⌋) n -⌊rn⌋ ⌊rn⌋ + 1
= β ⩽ 1, (b) α = β > 1, (c) α < β < 1; α close to β, (d) 1 > α > β; α close to β, (e) α ⩽ 1 ∧ α < β, (f) β ⩽ 1 ∧ α > β, (g) 1 < α < β, (h) α > β > 1.
π k2 π k1 = p 1 -p 2⌊rn⌋-n+1 . ( 39
)
When n is large, we obtain

π k2 π k1 ≈ p 1 -p n(2r-1)+1 . ( 40
)
Observe that the peaks have equal heights when p = 0.5, and when r = 0.5, the ratio is equal to p/(1p).

Besides, we have solved numerically by SCILAB the system of equations (35), which is possible for reasonable values of n, and obtained its solution π(k), k = 0, . . . , n. Table I shows the value of the ratio of the two peaks for various values of p, r as given by Eq. (39), compared to the output of SCILAB. Fig. 5 shows the computed distribution π for n = 100 for the one threshold case and also the case of the Beta distribution, compared to the histograms obtained from Monte Carlo simulations.

We comment on these results. The Markov approach permits to obtain the stationary probability distribution of the different states, for any value of n, without approximation. It is found that in the long run, even if any state has a nonzero probability to be reached, some as well as the values π k1 and π k2 computed numerically and the ratio between them (three rightmost columns).

states have a much higher probability than the others to appear. In the case of one threshold, we have analytically proved the presence of one or two peaks, and their positions when n is large perfectly coincides with what was predicted by the mean field approach. It is complementary to the results given by the mean-field approach, since the Markov approach is able to give the probability of occurrence of each stationary state. On the other hand, the complexity of the system of linear equations (35) induced by the Markov chain makes this approach not always tractable (e.g., with the Beta distribution). Nevertheless, we have shown that for reasonably large values of n (e.g., n = 100), this linear system can be solved numerically, giving a perfect fit with theoretical values, as shown by

VII. SUMMARY AND RESEARCH DIRECTIONS FOR THE FUTURE

In this paper, we investigated the threshold model with anticonformity under asynchronous update mode, which mimics continuous time. We considered two cases: (1) homogeneous, in which all agents had the same threshold and (2) heterogeneous, in which the thresholds are given by the beta distribution function. The homogeneous case with r = 0.5 is identical to the homogeneous symmetrical threshold model with anticonformity [START_REF] Nowak | Homogeneous symmetrical threshold model with nonconformity: Independence versus anticonformity[END_REF]. Moreover, it is almost identical to the majority-vote process [START_REF] Liggett | Interacting Particle Systems[END_REF][START_REF] Oliveira | Isotropic majority-vote model on a square lattice[END_REF]. The only difference between the models is when the number of active and inactive agents in the neighborhood of a chosen agent is equal. In such a case, the state of the system does not change within the majority-vote model, whereas within the threshold model the change is possible. From this point of view, the threshold model with anticonformity under asynchronous updating can be treated as a generalization of a majority-vote model.

On the complete graph, the homogeneous threshold model does not give particularly interesting results. The relationship between the stationary ratio of active agents and the probability of anticonformity consists of linear dependencies, similarly as for the homogeneous symmetrical threshold model [START_REF] Nowak | Homogeneous symmetrical threshold model with nonconformity: Independence versus anticonformity[END_REF][START_REF] Nowak | Promoting discontinuous phase transitions by the quenched disorder within the multistate q-voter model[END_REF]. The only interesting feature of this model is the discontinuity that appears at c = r = 1p for r > 0.5 and at c = r = p for r < 0.5. In the result, the system reaches one of two different steady states, depending on the initial conditions. Much richer behavior is observed in the heterogeneous model with thresholds given by the beta distribution function, parametrized by α, β, which allows tuning the model to the homogeneous one (α, β → ∞) and to maximally heterogeneous (i.e. described by the uniform distribution function).

A particularly interesting behavior is obtained if at least one of the shape parameters α or β is larger than one and both parameters are finite. In this case PDF has a shape that reminds those of actual trait manifestation in behavior, as reported by psychologists [START_REF] Fleeson | The implications of Big Five standing for the distribution of trait manifestation in behavior: Fifteen experience-sampling studies and a meta-analysis[END_REF], i.e., uni modal, not necessarily symmetrical, function with maximum at the value 0 < r < 1. In such a case a phase transition appears, which is continuous for α = β, and discontinuous otherwise. In the latter case, the transition involves phenomena typical of social systems, such as social hysteresis [START_REF] Scheffer | Slow response of societies to new problems: Causes and costs[END_REF] and the critical mass [START_REF] Centola | Experimental evidence for tipping points in social convention[END_REF][START_REF] Van Slyke | Perceived critical mass and the adoption of a communication technology[END_REF][START_REF] Mahler | The diffusion of interactive communication innovations and the critical mass: The adoption of telecommunications services by German banks Telecommunications Policy[END_REF].

The future research on the model can be conducted in several directions, related to the following questions:

• How the results would change if the threshold for anticonformity would be different than that for conformity? This question is inspired by the work on the q-voter model with generalized anticonformity [START_REF] Abramiuk-Szurlej | Discontinuous phase transitions in the q-voter model with generalized anticonformity on random graphs[END_REF]. In the q-voter model such a generalization resulted in switching from continuous to discontinuous phase transitions for some values of parameters. The question is if the same phenomena would be observed for the threshold model.

• How the structure of a network would influence the results? This question is inspired by the work on the symmetrical threshold [START_REF] Nowak | Promoting discontinuous phase transitions by the quenched disorder within the multistate q-voter model[END_REF]. It was shown that on random graphs with the degree observed empirically for social networks, the largest social hysteresis is observed for r ∈ (0.65, 0.85). This was a meaningful result from the social point of view and thus it would be desirable to check if it appears also in the asymmetric model studied here.

• How the results would change if the quenched approach to anticonformity would be used. In this version of the model, we used the annealed approach, in the sense that each agent could anticonform (with probability p) or conform (with probability 1p). However, we could use also the quenched approach, in which a fraction p of agents are permanently anticonformists. This question is inspired by the work on the q-voter model with nonconformity under quenched and annealed approaches [START_REF] Jędrzejewski | Person-situation debate revisited: Phase transitions with quenched and annealed disorders[END_REF]. It was shown that on the complete graph both approaches give the same result for the q-voter model with anticonformity, whereas different for the model with independence. The question is to what extend this result is universal.

FIG. 1 .

 1 FIG.1. Dependency between the stationary concentration of active agents c and the probability of anticonformity p for model with one threshold for different values of the parameter r (indicated in the plots). Solid lines represent stable fixed points obtained analytically from Eq. (16). Symbols represent Monte Carlo simulations from two initial conditions indicated in the legend. Arrows indicate the flow in the system. The length of each arrow correspond to the value of velocity dc

FIG. 4 .

 4 FIG. 4. Representative dependencies between the stationary concentration of spins up and the probability of anticonformity for model with beta distribution for different values of the parameters α and β:(a) α = β ⩽ 1, (b) α = β > 1, (c) α < β < 1; α close to β, (d) 1 > α > β; α close to β, (e) α ⩽ 1 ∧ α < β, (f) β ⩽ 1 ∧ α > β, (g) 1 < α < β, (h) α > β > 1.Solid and dotted lines represent stable and unstable steady states respectively, obtained with Eq. (27). The exact values of parameters in the plots are as follows: (a) α = β = 0.9, (b) α = β = 4, (c) α = 0.1; β = 0.2, (d) α = 0.2; β = 0.1, (e) α = 1; β = 3, (f) α = 3; β = 1, (g) α = 5; β = 8, (h) α = 8; β = 5. Symbols represent Monte Carlo simulations from two initial conditions, denoted in the legend. The results are averaged over 10 runs and collected after 10 4 MCS for system of size 10 4 .

  FIG. 4. Representative dependencies between the stationary concentration of spins up and the probability of anticonformity for model with beta distribution for different values of the parameters α and β:(a) α = β ⩽ 1, (b) α = β > 1, (c) α < β < 1; α close to β, (d) 1 > α > β; α close to β, (e) α ⩽ 1 ∧ α < β, (f) β ⩽ 1 ∧ α > β, (g) 1 < α < β, (h) α > β > 1.Solid and dotted lines represent stable and unstable steady states respectively, obtained with Eq. (27). The exact values of parameters in the plots are as follows: (a) α = β = 0.9, (b) α = β = 4, (c) α = 0.1; β = 0.2, (d) α = 0.2; β = 0.1, (e) α = 1; β = 3, (f) α = 3; β = 1, (g) α = 5; β = 8, (h) α = 8; β = 5. Symbols represent Monte Carlo simulations from two initial conditions, denoted in the legend. The results are averaged over 10 runs and collected after 10 4 MCS for system of size 10 4 .

FIG. 5 .

 5 FIG. 5. Stationary distributions of visited states for model with one threshold (upper row) with parameters (a) r = 0.5; p = 0.3, (b) r = 0.6; p = 0.7 and model with beta distribution (bottom row) with parameters (c) α = 8; β = 5; p = 0.1, (d) α = 3; β = 1; p = 0.4. Solid red lines are distributions obtained with Markov approach, black histograms are obtained with trajectories from Monte Carlo simulations for system of size n = 100 and thermalization time t = 2 • 10 6 MCS from 100 initial conditions evenly distributed on [0, 1] interval, averaged over 1000 independent runs.

TABLE I .

 I Example of results for n = 10 for different values of p and r under the condition p ⩽ r ⩽ 1-p. In the table are presented the theoretical ratio given by Eq. (39) (left column),

	p	r π k2 /π k1	π k1	π k2	ratio
	0.21 0.41	3.7619048	0.064177 0.241429 3.7619052
	0.21 0.61	0.0187835	0.2966234 0.005572 0.0187834
	0.25 0.5	0.3333333	0.218683 0.0728942 0.3333334

  Table I and with Monte Carlo simu-lations as well, see Fig 5.

ACKNOWLEDGMENTS

This research was supported by the National Science Center (NCN, Poland) grant number 2019/35/B/HS6/02530.