Jérémy Patrix
email: jeremy.patrix@unicaen.fr

Abdel-Illah Mouaddib
email: abdel-illah.mouaddib@unicaen.fr

Simon Le Gloannec
email: simon.legloannec@cassidian.com

Dafni Stampouli
email: marc.contatdafni.stampouli@cassidian.com

Marc Contat

Discrete Relative States to Learn and Recognize Goals-based Behaviours of Group

Keywords: I.2 [Artificial Intelligence]: Distributed Artificial Intelligence, I.5 [Pattern Recognition]: Models, Design Methodology, I.6 [Simulation and Modeling] Algorithms, Measurement, Security, Human Factors Agent-based simulation:Emergent behavior, Systems and Organisation: Complex systems, Self-organisation

come   L'archive ouverte pluridisciplinaire

INTRODUCTION

In this paper, we present a part of a system which has been developed within the EUSAS (European Urban Simulation for Asymmetric Scenarios) project. The project is financed by twenty nations under the Joint Investment Program Force Protection of the European Defense Agency (EDA). Within this project, we aim to obtain improved situation awareness by using new methods for detecting high-level information from observed variables of the monitored crowd. Previous work [START_REF] Patrix | Detection of primitive collective behaviours in a crowd panic simulation based on multi-agent approach[END_REF] developed methods to assess risks during social gatherings in the case of emergency and panicked crowd evacuation, by identifying problematic escape paths and possible stampedes and points of crowd crushing. The first methods allow training, recognition and anticipation, in real-time, of the recurrent emerging collective movements of panicked agents (each agent corresponds to a person in the crowd). Then, in order to detect the reasons of this type of behaviors, we developed methods for group detection. This paper presents consequent work on the detection of these groups' behavior, their goals and intentions, and the detection of the reasons that caused these behaviors. Group behavior modeling methods were used, such as Hidden Markov Model (HMM) by learning from observations of successive states of members and which extract high-level facts of current behaviors as Complex Event Processing (CEP).

Real-time detection of a multi-agent behavior is a highly complex problem. Indeed, if the group state is defined based on the states of its members, a model of group behavior will have a number of the variables |A|.|S|.t where A is the set of agents and S is the number of variables for each agent state from the start of the observation till the current time step t. Therefore, the size of the behavior model becomes intractable, over time.

Our approach consists in defining a set of discrete relative (DR) state sequences that model goal-based group behavior. The set of DR states represents a discretization of the evolution between two successive times concerning the difference between a goal state and a centroid state of a group. This centroid of the group is computed from local states of people of the group. By this way, the number of agents/people in the group has no impact on the complexity of the recognition problem of a collective behavior.

An important fact about the goals of a group is that they are all entities of the environment (excluding the members of the current group, such as other people, other groups, locations, critical infrastructures). The action sequence of a group has been decided according to these external goals/entities, for example: a group is following this target or a group is fleeing from an identified threat. The proposed method provides a probability value denoting the likelihood that this group is executing this behavior according to the specific goal, thus they are called goal-based group behaviors.

The algorithm is used to obtain the most likely groups and goals among observed agents. The probability that such group executes such collective behavior to reach a specific goal is established according to an observation sequence and a goal-based group behavior model. During the simulation, we used a highly realistic virtual environment, as presented in the last section.

Section 2 and 3 present the background and previous work relevant to our framework, followed by our contribution in section 4. Finally, section 5 presents the results from the experiments in the virtual environment.

BACKGROUND

In order to facilitate the comprehension of our model and the suitability of our experiments, this section presents relevant background knowledge in the areas of teamwork analysis.

The objective of the teamwork analysis [START_REF] Sukthankar | Activity recognition for agent teams[END_REF] is to identify the team's goal, the executed actions of members to achieve the goal and the coordination patterns between its members. Three general classes of characteristics are then potentially valuable for the team's activity/plan recognition [START_REF] Sukthankar | Activity recognition for agent teams[END_REF]: (1) spatial relationships between team members and/or physical landmarks that remain fixed over a period of time; [START_REF] Bowman | The dual node network (dnn) data fusion & resource management (df&rm) architecture[END_REF] temporal dependencies between behaviors in a plan or between actions in a team behavior; and (3) the coordination constraints between agents and the actions that are performed. Some properties of the observed interactions are invariant over all possible permutations of behavior, but the multiagent activity/plan recognition includes two difficulties [START_REF] Sukthankar | Activity recognition for agent teams[END_REF] : (1) the multi-agent behaviors are more variable in the observed spatiotemporal traces than simple behavior; (2) the team composition can change over time.

The observation sequence of agent behavior can be segmented to obtain its state sequence. In real-time, this sequence is continuous and thus must be tested at regular intervals. This requires additional memory for the state space of observed data for the learning and additional execution time for the methods that learn and recognize the behavior models.

In our model, we use Hidden Markov Models (HMM) that allows representing a behavior on a set of hidden states with a transition probability distribution between these states. All observable state sequences of this behavior allows us to obtain the transition probabilities.

We adapt the HMM definition [START_REF] Rabiner | An introduction to hidden markov models[END_REF] by this 5-tuple Q, A, O, B, Π :

Q is a finite set of hidden states, where qi is the i th hidden state;

A is a finite set of transition probabilities, where aij = P r(st+1 = qj|st = qi) is the transition probability that the hidden state will be qj after qi ;

O is a finite set of observable states, where oi is the partially observable state of qi (as obtained from a monitored system after using data fusion methods to reduce the noise);

B is a finite set of observation probabilities, where b ik = P r(o k |qi) is the probability to observe o k in the state qi;

Π is a finite set of initial state probabilities, where πi = P r(s0 = qi) is the probability to be initially in the state qi. and allows three main computations [START_REF] Rabiner | An introduction to hidden markov models[END_REF]:

1. The probability of the current state by the equation P r(st = qt|Ot) that gives the most likely hidden state qt at the time t knowing an observation sequence Ot = {o1, ..., ot}. It is solved by a dynamic programming algorithm such as Forward-Backward [START_REF] Rabiner | A tutorial on hidden markov models and selected applications in speech recognition[END_REF];

2. The most likely hidden state sequence knowing an observation sequence. It is solved by the Viterbi algorithm [START_REF] Rabiner | A tutorial on hidden markov models and selected applications in speech recognition[END_REF];

3. The HMM learning is able to give the maximum probability when it estimates a sequence. It is solved by the Baum-Welch algorithm [START_REF] Rabiner | A tutorial on hidden markov models and selected applications in speech recognition[END_REF] that uses the Expectation-Maximization method.

In our work, different levels of observations and modeling produce a stochastic representation of a group behavior: (1) the micro-level is about the roles of the agents impacting on their possible individual activities; (2) the meso-level is about the activities executed in teamwork; (3) the macro-level is about the goals that explain the interactions between each group and the environment.

To model a behavior of agents, the model must be canonical [START_REF] Sukthankar | Automatic recognition of human team behaviors[END_REF] whatever the position, the speed, the orientation, the observation scale, the number of agents and goals. For a canonical representation of the team state, the concept of centroid is used [START_REF] Sukthankar | Automatic recognition of human team behaviors[END_REF]. This corresponds to a vector where each variable is the average value of states variables of the team members. Thus, our team's behaviors are defined on their relative movement and not on the absolute variables of agents to be recognized in any monitored situations. The extracted characteristics of the behavior models could be classified by minimal information and discrete according to their type of agent, environment or team. Such examples are the Centroid-Relative Position Vector [START_REF] Luotsinen | Role-based teamwork activity recognition in observations of embodied agent actions[END_REF] or the Role-Relative Position Vector that change according to members' roles during team behavior.

p c x,y,t = |c| i=1 p a i x,y,t |c| (1)
At time t, the group position p c x,y,t of the centroid state s c t is computed by the average of positions p a i x,y,t of each agent a i of the group c. The centroid state s c t includes all required information of the agents (such as position and speed).

In conclusion, a centroid state formalizes a group state with an unspecified number of members. The sequence of these states characterizes team behavior, and the behavior can be learned and recognized using HMM methods, thus achieving detection of group behavior. Our model presents the innovative points that consist of a discrete relative state sequences which are used to formalize and detect goal-based group behaviors. Before the proposed method is explained in greater detail, previous related work is presented in the following section.

RELATED WORK

There are several studies that investigated the requirements, advantages and limitations of methods on group behavior and goal detection. This section describes the evolution of relevant technical solutions to demonstrate the advantages achieved by the proposed model.

The behavior planning and detection have a difference in the direction of information propagation like in the Dynamic Bayesian Networks [START_REF] Giersich | Towards an integrated approach for task modeling and human behavior recognition[END_REF] that propagates it in the two directions and model time influence. The probabilistic plan/activity recognition [START_REF] Robert | A bayesian model of plan recognition[END_REF][START_REF] Shi | Propagation networks for recognition of partially ordered sequential action[END_REF] of Bayesian Network (representing links between events/actions with a partial order constraint) is applied by an inference algorithm exactly as the Viterbi algorithm on HMM (computing the posterior probabilities of distribution on possible interpretations).

In order to model complex behaviors, the main idea has been the decomposition of the parameter space to decrease their complexity: in initial, accept, intermediate and reject states (Behavior HMM [START_REF] Han | Automated robot behavior recognition[END_REF]); or until each observed agent characteristic (Observation Decomposed HMM [START_REF] Liu | Multi-agent activity recognition using observation decomposed hidden markov model[END_REF]); or for each distinct temporal granularity level and each kind of information (Layered HMM [START_REF] Oliver | Layered representations for human activity recognition[END_REF]); or (Switching Hidden Semi-Markov Model [START_REF] Duong | Activity recognition and abnormality detection with the switching hidden semi-markov model[END_REF]) using a duration distribution for each state to consider that a state remains unchanged for a certain duration before its state transition. Group behavior can be executed by a variable number of agents, thus the decomposition solution increases the complexity of the problem (because each group behavior should be decomposed and learned for each possible team size and roles members).

In order to detect team behaviors, the probabilistic recognition of agent's role can be executed before the team activity recognition to improve the detection quality. The Unique or Multiple Role Assignment algorithm [START_REF] Luotsinen | Role-based teamwork activity recognition in observations of embodied agent actions[END_REF] selects agents for each identified valid roles having the highest probabilities (travelling from the root to the leafs of roles in a probabilistic decision tree). For the team activity recognition, they represent an Idealized Team Action [START_REF] Luotsinen | Role-based teamwork activity recognition in observations of embodied agent actions[END_REF] after the normalization of observations to deal with the noise. This method's limitation is that it must produce a library of the same team activity for each member condition, as in the Multi-Agent Plan Recognition [START_REF] Banerjee | Multi-agent plan recognition: Formalization and algorithms[END_REF]. In order to detect a team's behavior, this method identifies the structures of agent behaviors in the team by a matrix, where pij is the expected action of the member of the line j at the step of the column i and must be compare in each team plan P .

Recent works have been focused on Abstract HMM [START_REF] Bui | Policy recognition in the abstract hidden markov model[END_REF] models which are used to recognize an agent policy on long term. Each upper level abstract behavior includes lower level simple action sequences. Therefore to manage its explosion in the state space, a hybrid inference method of Factored State-AHMM [START_REF] Tran | Factored state-abstract hidden markov models for activity recognition using pervasive multi-modal sensors[END_REF] is used to remove the redundant probabilistic links of state transitions.

Each state of each level of a Hierachical HMM [START_REF] Bui | Hierarchical hidden markov models with general state hierarchy[END_REF] can be itself a HHMM. This hierarchical policy recognition of a team is based on a multi-agent decision method called Hierarchical Multiagent Markov Process [START_REF] Saria | Probabilistic plan recognition in multiagent systems[END_REF] that coordinates actions of agents on a upper level by a controller and without coordination on a lower level. An end state for each level indicates when the activity is finished in order to go to the upper level state.

Using the Joint Probabilistic Data Association Filters method, HHMM-JPDAF [START_REF] Nguyen | Recognising behaviours of multiple people with hierarchical probabilistic model and statistical data association[END_REF] recognizes the behaviors of multiple persons on its set of HHMM and the noisy observations. Its learning algorithm discovers [START_REF] Nguyen | Discovery of activity structures using the hierarchical hidden markov model[END_REF] the HHMM structure and estimates the link probabilities by the Expectation-Maximization method.

For recognition on long term, a memory node of each level in the Abstract Hidden Markov mEmory Model (AHMEM) [START_REF] Bui | Efficient approximate inference for online probabilistic plan recognition[END_REF] memorizes the actions sequence that is executed according to the chosen activity by the current policy of the upper level. Knowing a noisy observation sequence, an approximate probabilistic inference computes, in real-time, the probability of various possible behaviors of AHMEM [START_REF] Nguyen | Learning people movement model from multiple cameras for behaviour recognition[END_REF].

In addition to a real-time probabilistic behavior recognition of observed multi-agent system, M(ulti-agent)-AHMEM [START_REF] Gaitanis | Cooperative team behaviour recognition for multimodal fusion[END_REF] allows a contextual, non-deterministic approach and a high level multimodal fusion of various sources using an ontology as library of behaviors describing the possible relationships between entities and the environment.

The learning of these last multi-agent behavior models is more and more complex according to the multiple kinds of states. Furthermore, the use of approximate inference methods for a real-time detection decreases the quality of detection. This raises the question whether there are alternative ways to HMM.

Techniques based on learning could be an alternative to HMM. Hierarchical Conditional Random Field [START_REF] Liao | Hierarchical conditional random fields for gps-based activity recognition[END_REF][START_REF] Vail | Conditional random fields for activity recognition[END_REF] encodes relations between observations including directly the conditional distribution, or Inverse Reinforcement Learning (IRL) [START_REF] Choi | Inverse reinforcement learning in partially observable environments[END_REF] describes the intention/goal/behavior of each agent by computing a reward function to reproduce the observed policy of the agent. However, these techniques suffer from the need of significant data/time for learning and are not suitable to detect goal-based group behaviors.

Like an alternative to the Bayesian inference, the Dempster-Shafer theory is useful to assign any agent to a team and behavior for each time step. The STABR (Simultaneous Team Assignment and Behavior Recognition) algorithm [START_REF] Sukthankar | Simultaneous team assignment and behavior recognition from spatio-temporal agent traces[END_REF] creates a set of potential teams. The formation recognition by the RANSAC (Random Sampling and Consensus) algorithm generates assumptions of teams that explain all the observed space-time traces of members. The advantages are the quick behavior recognition in comparison with HMMs and the team member reassignment. But it cannot guaranty obtaining the best match between the observed data of the team behavior and the selected behavior models among all assumptions.

Another possible solution is the Complex Event Processing (CEP) [START_REF] Luckham | The power of events: an introduction to complex event processing in distributed enterprise systems[END_REF] that includes an inference engine that uses the forwardand/or backward-chaining on its rules to detect the behaviors and situations requiring a reaction. CEP can be represented using levels of the JDL Data Fusion model [START_REF] Bowman | The dual node network (dnn) data fusion & resource management (df&rm) architecture[END_REF]: (1) the standardization rules on various sources, (2) the filter rules and the simple rules execution on Complex Events (CEs), (3) the aggregation rules for new CEs and the detection rules for significant items in an events cloud, (4) a feedback loop. But the inference engine has difficulties to deal with a great volume of facts with regular updates by many CEs and production rules. We have already implemented and often used this technique and it will be used as the comparison method with our model during experimentation.

In conclusion, three main approaches are presented in this state of the art section: probabilistic, grammatical and statistical / fuzzy methods. All these approaches are unable to detect in the same way the collective behavior and the target goal. The following section presents the proposed approach based on an aggregated state (centroid) of the group allowing the recognition of the behavior and the goal of the group, in real-time.

THE MODEL OF GOAL-BASED GROUP BEHAVIOR

In the framework of the EUSAS project, an algorithm, similar to the one presented in [START_REF] Chang | Probabilistic group-level motion analysis and scenario recognition[END_REF], detects groups of agents. This study advances the previous work in order to achieve goal-based group behavior detection. In this section, we present how to compute a discrete relative (DR) state, how to model a collective behavior according to a goal and how to detect goal-based group behavior. In the study, the state of an agent contains only its position p and its speed v.

Discrete Relative states

To compute a DR state as in Figure 2, we get observations {o a i t , o g j t } on the hidden states {s a i t , s g j t } for each agent a i of the group c and each goal g j / ∈ c (which could be mobile) from the time 0 to t. Each ot and st contains the position and the speed pt, vt .

1. We compute {s c t-1 , s c t } the centroid states at time t -1 and t of a group c using equation 1.

2. We compute relative states {∆s c,g t-1 , ∆s c,g t } on the difference between centroid states {s c t-1 , s c t } and goal states {s g t-1 , s g t }. The relative variable ∆r c,g t ∈ ∆s c,g t is given by:

∆r c,g t = |r c x,t -r g x,t | 2 + |r c y,t -r g y,t | 2 (2)
where rx,y,t is a spatial variable (as px,y,t of the equation 1).

3. We discretize the relative states {∆s c,g t-1 , ∆s c,g t } on their difference from t -1 and t for each variable ∆r c,g ∈ ∆s c,g . The discrete relative variables {drn, dr d , drc, dri} are introduced to represent any DR state which are, in our case, the DR positions {pn, p d , pc, pi} and the DR speeds {vn, v d , vc, vi}. These DR variables are given by: drn (is null) ∆r c,g t < δr min dr d (decrease) if ∆r c,g t < ∆r c,g t-1 -δrmin dri (increase) ∆r c,g t > ∆r c,g t-1 + δrmin drc (is constant) |∆r c,g t -∆r c,g t-1 | < δrmin The value δrmin indicates a minimal state transition between two successive times according to the variable r ∈ s, thus it is directly proportional with the number of obtained observations per seconds. If we take the relative position p and speed v between a group and a goal from t -1 to t, the computed set of DR states represents:

Q DR =     pn,
   
The advantage of this discretized state space is to facilitate the HMM learning and this kind of discretized state is human-readable.

The DR state pn, vn represents a null discrete relative position pn and a null discrete relative speed vn (the group position and speed is equivalent to the goal position and speed) from t -1 to t.

In a sequence of DR states, we can show the relative movement over time between a group and a goal. For example, the three following DR state sequences show three different goal-based group behaviors:

∆S c,g 0,4 = { pn, vn , pn, vi , pi, vi , pi, vc , pi, vc } ∆S

Goal-based Group Behavior Assignment

To produce goal-based group behaviors detection, we consider the tuple O c,G t , Λ, F where:

O c,G t = {{o a 1 0 , o a 2 0 , ..., o a |c| 0 , o g 1 0 , ..., o g |G| 0 }, ..., {o a 1 t , ..., o g |G| t
}} is an observation sequence of each agent of the group c and each goal of G from 0 to t, where {o a i t , o g j t } is an observation of the i th agent and the j th possible goal that an observed group c can seek to reach; Λ is a set of HMMs, where λ h ∈ Λ is the HMM modeling the h th defined group behavior;

F is a set of functions, where f h is a function (i.e. section 4.1) that computes the DR state sequence ∆S h t from an observation sequence according to HMM λ h . P r(λ|Ot) is the probability that the behavior λ is executed during the observation sequence Ot (computed by Forward-Backward algorithm [START_REF] Rabiner | A tutorial on hidden markov models and selected applications in speech recognition[END_REF]). During the T last observations (in order to capture the recent behavior), the computation of observations sequence O c,g j t-T,t in DR states sequence ∆S c,g j t-T,t is given by:

M c t (λ h , g j ) = P r(λ h |∆S c,g j t-T,t = f h (O c,g j t-T,t )) (3)
For each group behavior λ h and each possible goal g j , we compute the probability (3) that a group c is executing λ h to reach g j . The result gives us thus the following matrix:

M c t =        M (λ 1 , g 1 ) • • • M (λ 1 , g j ) • • • M (λ 1 , g |G| ) . . . . . . . . . . . . . . . M (λ h , g 1 ) • • • M (λ h , g j ) • • • M (λ h , g |G| ) . . . . . . . . . . . . . . . M (λ |Λ| , g 1 ) • • • M (λ |Λ| , g j ) • • • M (λ |Λ| , g |G| )       
We obtain the most likely tuple λ H , g L by:

M (c, t) = argmax ∀g j ∈G ∀λ h ∈Λ {M c t (λ h , g j )} (4)
Here, only one goal-based group behavior is selected for the group by the best probability of M c t . But a group can execute different group behaviors for different goals at the same time. We obtain the most likely group behavior λ H for each goal ∀g j ∈ G by the result of:

M (c, g j , t) = argmax ∀λ h ∈Λ {M c t (λ h , g j )} (5)
In the same way, for another use case, for a group, we could obtain M (c, λ i , t) the most likely goal g h for each group behavior λ i on M c t . To obtain a probability on the recent behavior of agents, we keep only the T last observations. We will show in experiments the impact on the detection quality by changing the observation sequence length T for the learning and the recognition.

In order to understand, at this point, following the instructions of this section, we model a goal-based group behavior using discrete relative state sequences computed from observation sequence of monitored groups of agents and their possible goals. The following section presents our method to extract the specific goals that are responsible of observed behavior sequences of groups.

Valued Goal-based Group Behavior

In a monitored situation with lots of possible goals for a group, we detect its current goal-based group (GbG) behaviors. In order to detect if the group is executing a sequence of GbG behaviors according to a particular goal (its intention), we measure its value using the following process that we call Valued Goal-based Group Behavior Detection.

This value can be measured in different ways according to the kind of intention that we want to detect. From the time t -T to t, we use the following goal values in order to represent:

• V c i ,g j t-T,t the intention of the group c i according to the goal g j ;

• V c i ,G j t-T,t the intention of the group c i according to a set of goals G j (∈ G, for example: all the security forces, the critical infrastructures, the safe areas or the dangerous areas...);

• V C i ,g j t-T,t the intention of a set of groups C i according to the goal g j .

In order to measure this goal value, we use a distribution where a goal g and a GbG behavior λ h of a group c give the value V (c, g, λ h ). A distribution is created according to the intention φ p (∈ Φ the set of defined intentions) that we want to detect by the observed GbG behavior sequence. The advantage of a distribution is that it can be created or learned to detect each type of intentions that is needed. So, the equation of the goal values are:

V φ p ,c i ,g j t-T,t = t tc=t-T V φ p (c i , g j , λ h = M (c i , g j , tc)) (6) V φ p ,c i ,G j t-T,t = g l ∈G j V φ p ,c i ,g l t-T,t ; V φ p ,C i ,g j t-T,t = c l ∈C i V φ p ,c l ,g j t-T,t
For example, in order to detect the civilians with a dangerous intention φ p for themselves: the distribution gives a positive value when civilians move towards a safe area and a negative value when the groups move toward a dangerous area, and after all groups executing φ p will have a higher positive goal value. According to the distribution of the intention φ p ∈ Φ, the group c is moving toward the specific long-term goal g φ p ,c is:

g φ p ,c t-T,t = argmax ∀g j ∈G V φ p ,c,g j t-T,t (7)
According to the monitored situation and what we search to identify, a valued goal-based group behavior indicates the intention of the group on the long-term on two levels of observation: firstly, a general information about its specific goal, and secondly, if its intention is normal (friend) or abnormal (foe). Thus the ability to learn and recognize valued goal-based group behaviors using HMM methods in real-time is shown during the following experimentation section.

EXPERIMENTS

During the experiments, the proposed method of the relative states in HMM is compared with CEP. Our implementation in java communicates with Virtual Battle Space 2 (VBS2)1 a highly realistic simulated environment and can manage HMMs and CEPs. For the HMMs, it uses a jahmm 2 library available on-line (managing states, observations and standard algorithms relating to HMMs). It uses also Drools a CEP rule engine (leading java open source rule engine federated in Jboss). During a scenario in this software, it extracts essential information of agents to apply the proposed model of detection and it displays the results in real-time.

This section begins with the HMM learning of behaviors. The following subsection presents the results from the experiments in order to show how the proposed model can be used to help in a crisis management situation. This section also includes an experiment of a simulated crowd during a social gathering with an human carrying a bomb that moves toward a press reporter (in order to detect the abnormal behavior of the crowd).

The HMM learning of behaviors

For the proposed use case, between a group and a goal, the probability of four different GbG behaviors needs to be determined in real-time. These behaviors are: λ c→g moving towards a goal, λ c←g moving away from a goal, λ c g moving around a goal and λ c g not moving according to a goal.

In order to define HMM on a goal-based group behavior, a learning model on the observation data resulting from various scenarios is created using the Baum-Welch algorithm [START_REF] Rabiner | An introduction to hidden markov models[END_REF]. This algorithm enables us to maximize probabilities of each given observation sequence by updating the parameters A, B and Π (i.e. section 2).

During the simulated scenarios, the teams are composed of 4 to 16 civilians that move in group formation. There is some (mobile or not) goals (defined by other agents or specific positions). The environment has also few obstacles that the teams avoid. The observations of each position and speed of agents is obtained each second to compute their DR state sequences (i.e. the process is described in section 4.1). . According to each behavior, only a subset of all relative states is used for learning and the transition probabilities are different even if the two behaviors overlap. So, we can learn various different GbG behaviors using our method.

HMM versus CEP

In order to compare the results of the proposed discrete relative states on Hidden Markov Model, the Complex Event Processing was used as an alternative method because it is often used for the complex event detection [START_REF] Bowman | The dual node network (dnn) data fusion & resource management (df&rm) architecture[END_REF][START_REF] Luckham | The power of events: an introduction to complex event processing in distributed enterprise systems[END_REF]. The following experimentations of detection are produced on a set of scenarios where collective movements are simulated by a multi-agent system in VBS2.

Drools was used to write rules that take observations of agents in low level facts and extract high level facts of collective behaviors and goals of groups. These rules compute the same function (as explained in section 4.1 and 4.2) to create high level facts including the same result of the centroids, the DR state sequences and M c t (the most likely collective behavior for each goal for a group).

The illustrations of Figure 4 are examples of monitored situations used to test the precision of our detection method. The dotted lines show the paths used by the teams (black for the centroid of each team and grey for each member). Each group has its own color and a red line links each member of the centroid. For each group during the detection, note that the other entities (alone agents, groups or important positions) compose the set of its possible goals and are presented by blue diamonds.

During the experimentations, the behavior detection quality is measured by the percentages of the precision expressed by A/B: where A is a value that increases one by one point when a detection corresponds correctly to the observation, and B is a value that increases one by one point when the system obtain a new observation that must be detected.

With the four GbG behaviors λ c→g , λ c g , λ c←g , λ c g , we added an inconsistent GbG behavior λ c g that is a HMM by default without learning. If λ c g has the highest probability then the observation sequences used during the behavior learning require further corrections. 1: The precisions of Goal-based Group Behavior Assignment on multiple scenarios (i.e. Figure 4) using CEP and HMM (according to the observation sequences lengths used to learn and recognize them).

sequence length of 8 as we found that for higher length the precision decreases. In order to obtain the best detection quality, the last five observations (one per second) in the tested sequences (of these scenarios) always give the highest precision. The decomposition of sequences of GbG behaviors to learn them with only two observations give a better precision. The solutions with a short length of observation sequences (as 2,2 in the table 1) detects more quickly each new GbG behavior but each new noisy observation has more impact (e.g., when the group avoids an obstacle during their observed movements). The solution with a long length of observation sequences (as 5,5) detects the current GbG behaviors even with more noisy observations but each new GbG behavior takes longer to be detected (after two new observations). The majority of errors during the detection occurs when a group changes its current GbG behavior. In order to increase the precision more than 92%, we added each possible GbG behavior (GbGB) that overlaps the evolution of an old GbGB toward a new GbGB. But by increasing the size of the library of GbGB models, we add a difficulty for an human user to understand what is happening.

CEP updates instantaneously the current state of the GbG behavior but without taking account the uncertainty of observations. In order to help in a crisis management situation, the poor precision (58%) indicates that the behaviors of groups and their goals can not be detected by CEP. On the opposite, the Goal-based Group Behavior Assignment reduces the impact of observations under uncertainty and always detects the nearest goal-based group behavior after two observations.

As a conclusion the precision of our proposed method can be used to help in crisis management context, in particular with the valued goal-based group behavior detection used as demonstrated in the following experiment.

Threat detection on complex dynamic environments

In the framework of the EUSAS project in VBS2, we simulated a busy urban street. Two hundred people can be seen walking peacefully in a road until they realise that there is a suicide bomber among them. They can see a man wearing a concealed bomb under his shirt. Then, the people in crowd start to flee in panic. However, there are people that do not see the man, and without realising the imminent danger they continue moving towards him. At the scene there is also a reporter that observes the situation and tries to understand what is happening. A screenshot of the simulated street is shown in Figure 5(a).

The proposed methods are applied to the simulated scenario. Initially, the crowd was separated into groups according to their observed state of the group. Thus we applied the valued goal-based group detection (i.e. section 4.3) in order to detect the intention φ ⇔g j ⇒ of groups run away from a specific entity (that is a nondetectable threat) within the crowd and the intention φ c i -→g j of a group than continue its movement along a long-term predefined goal. The result of the group detection is shown in Figure 5(b).

In order to identify φ ⇔g j ⇒ , we verify for each goal, if its goal value V φ ⇔g j ⇒ ,C i t-T,t

is not abnormaly high using a distribution that gives a positive value when groups move around or away from the goal, else it gives a negative value: V ({λ

C i ←g j , λ C i g j , λ C i g j , λ C i g j , λ C i →g j }) = {1, 1 2 , 0, - 1
2 , -1}. And φ c i -→g j the longterm goal g j of the group c i is identified by the highest goal value

V φ c i -→g j t-T,t using the distribution: V ({λ c i →g j , λ c i g j , λ c i g j , λ c i g j , λ c i ←g j }) = {1, 1 2 , 0, 1 2 , -1}
. With these each non-searched goal value is near to zero over time, and the searched goal value is quickly higher than the T value (because of is computed from t-T to T ) and increased by the number of goals or groups (when the goal value is the sum of their goal values, i.e. Equation 6).

Figure 5(a) shows the panicked groups running away for the threat. Note that not everyone realised the threat and some people continue along their original path. Figure 5(c) shows the results of the detections. It presents an overview of the situation for the analyst. First the threat is detected based on the behavior of the people towards this entity. Also the groups that did not detect the threat and are in danger can also be detected. This enables the security forces to predict the groups in danger and from the group members predict the number of casualties.

Another important aspect of the method is that if the suicide bomber tries to hide within the group the system treats the entire group as dangerous, rather than allowing the threat to be hidden within the majority of an innocent group. Therefore the system can continue tracking the danger even if he tries to conceal himself by merging with another group. However, a difficulty that the method faces is to create a correct distribution of goals values according to the library of GbG behaviors because each abnormal intention executes specific behaviors that must be valued correctly in order to avoid fake detections. This provides a powerful tool to the security forces to understand the situation and predict threat in order to take better informed decisions on the developing situation.

This could involve resource allocation to investigate further the detected threat, or eliminating the threat once it is confirmed. Furthermore, it could provide information about the member of public in greater danger so that the security forces can focus on accessing them and protecting them. Finally, in case there are critical infrastructures and vulnerable entities, the system can be used to assess (a) People that realised the presence of the suicide bomber can be seen running and fleeing the scene. The people that did not see the bomber can be seen continuing their movement or looking surprised at the panicked people.

(b) A 2D view of the situation. The marks denote the different group detected and their direction of movement.

(c) It is easier to follow each centroid movement of groups. The valued goal-based group detection reveals that there are groups that avoid and move away from a single agent during its movement toward the player. the intention of the target towards these assets, based on whether he is moving towards them or not.

CONCLUSION

We have shown that the use of discrete relative states in Hidden Markov Model may successfully detect, in real-time both multiple collective behaviors and goals of groups based on observations. We have compared the proposed approach with a different method called Complex Event Processing, and found the results are superior in terms of precision and tolerance to noises.

In future work on asymmetric threat environments, we plan to apply our approach on new abnormal collective behaviors such as panic and violent crowd behavior by the fusion of new types of information. The automatic identification and learning of new collec-tive behaviors including agent roles is envisioned to obtain a fully automatic collective behavior detection system. This last objective could allow an analyst to make an informed decision in the case of an asymmetric threat, taking into account the detected results.

Figure 1 :

 1 Figure 1: An example of HMM with two hidden states and three observable states This system is kept under following constraints ∀i ∈ Q,

c,g 5 , 9 =

 59 { pc, v d , pc, vi , pc, vc , pc, vc , pc, v d } ∆S c,g 10,14 = { p d , vc , p d , vc , p d , v d , pn, v d , pn, vn } First during ∆S c,g 0,4 , the group c quitted the goal g and has moved away from it, second during ∆S c,g 5,9 , c has moved around g at the same distance, and last during ∆S c,g 8,12 , c has moved towards g and is now close to it.

Figure 2 :

 2 Figure 2: A discrete relative state computed from a set of observable states of n agents of a group c and a goal g An observed sequence of this DR states can be seen as a low level group behavior sequence of a recurrent high level group behavior. Each group's behavior has its own set of recurrent state sequences and we use HMM to discover it.

  (a) The team is moving towards the goal. (b) The team is moving away from the goal.

Figure 3 :

 3 Figure 3: Two HMM examples of goal-based group behaviors that have been learned (and they have their own transition probabilities between the discrete relative states). Two learning results are shown on Figure 3. Each of 16 DR states (based on the relative positions and the relative speeds) has a

Figure 4 :

 4 Examples of monitored experimentations. transition probability not null to go to another DR state. Edges represent their transition probabilities (with a proportional size only if they have a probability Aij > 1/|Q DR | = 0.0625, ∀i, j ∈ Q DR ). Each double circle indicates us the possible initial states of the behavior. The GbG behavior of Figures 3(a) and 3(b) are completely opposed and we can see that their initial and final DR state are almost at the opposite (as pn, vn that is the final DR state of 3(a) and the initial state of 3(b))

Figure 5 :

 5 Figure 5: A situation with a high-level risk.

Table 1

 1 presents the tests for the GbG behavior learning and recognition on different lengths of observations sequences from 2 to ∞ (by an iterated method). The table presents results up to a 2% 84.0% 83.5% 83.5% 83.3% 83.6% 83.6% 3 86.2% 84.6% 84.2% 84.2% 83.7% 83.8% 83.8% 4 87.0% 85.4% 85.8% 86.2% 86.6% 85.3% 83.8% 5 88.1% 88.1% 85.8% 85.9% 85.9% 85.7% 86.1% 6 85.4% 85.1% 85.0% 85.0% 85.2% 82.7% 82.6% 7 81.8% 82.2% 81.3% 81.5% 82.2% 82.3% 82.3% 8 77.4% 78.7% 78.2% 78.3% 78.2% 78.6% 78.8%

	Length of recognized sequences 2 86.Table 1 CEP HMM 2	3	58,06 % Length of learned observation sequences 4 5 6 7	8

Bohemia InteractiveSimulations, Australia,

2011, http://vbs2.com 2 J.M. Francois, Jahmm (Java Hidden Markov Model), 2006, http://code.google.comp/jahmm/


Keywords: queue, Phase type distribution, Skewness, Approximate solution, Relative error

High variability of system parameters is a complicating factor in the modeling of the performance of big data systems. In this paper, we assess the potential inaccuracy of several existing approximations for evaluating the mean number of jobs queued in a parallelized device that can be represented as an M/Ph/m queue. Unlike existing studies, we consider the effect of the third moment of the service time, or equivalently, its skewness. We show that the approximations accuracy can be poor even for "easy" examples with a low coefficient of variation of the service time. Our examples demonstrate the important influence of the skewness of the service time distribution on the accuracy of the approximations. None of the approximations accounts for this property. We provide recommendations for the choice of the approximation that allow the user to choose the best suited approximation based on the actual queue parameters.

A note on the accuracy of several existing approximations for M/Ph/m queues Thomas Begin * and Alexandre Brandwajn †

I. INTRODUCTION

With the rapid emergence of cloud computing, the proliferation of social networking sites and the continuous growth of storage systems, practitioners are increasingly facing the challenge of dealing with extremely large datasets, commonly referred to as big data. The processing of these datasets raises many new issues related to their capture, storage, search, sharing, analysis, visualization and performance evaluation. For this latter purpose, the process of modeling systems, which include big data, is generally seen as complex due to, among other things, the high variability of some system parameters. For instance, in the case of a high-performance disk controller, the time needed to process an I/O request is likely to exhibit a skewed distribution because of the underlying multi-tier storage architecture.

Thus, when modeling a system, practitioners are often inclined to replace complex variables with more compact but approximate descriptions. By doing so, they tend to simplify the model and hopefully its solution. However, the resulting cost in terms of (in)accuracy is often unclear. In this paper, we assess the potential inaccuracy of several well-established approximations for evaluating the performance of a parallelized device. We focus our study on the case of highly variable workloads as it is likely to be the case in the context of big data.

The remainder of the paper is organized as follows. In Section II we present the simplest model to represent the behavior of a parallelized system. We briefly review existing approximations to evaluate its steady-state performance, and we detail our methodology to assess the applicability domain of these approximations. Section III presents the numerical results and discussions. Section IV concludes the paper.

II. MODEL OF PARALLELIZED DEVICE AND ITS SOLUTION

The simplest yet meaningful model for a parallelized device with a skewed service time distribution is the so-called M/Ph/m queue. This queue represents a system with m homogeneous servers, where the time between arrivals are exponentially distributed and job service times have a phase type distribution 1 . The buffer (i.e., queueing room) size is assumed to be infinite, so there is no limit on the number of jobs it can contain. The arrivals occur at a rate λ. The service time distribution has a mean of 1/µ and a coefficient of variation of c B . Recall that the coefficient of variation is defined as the ratio of the standard deviation to the mean, and can be seen as a normalized second order moment. The servers utilization is given by ρ = λ/(m.µ). Of course, if ρ < 1, i.e., λ < m.µ, the model is stable and the number of jobs in the system has a stationary distribution. Figure 1 depicts the M/Ph/m queue considered while Table I summarizes the principal notation used in our paper.

In the general case, no simple derivation exists to compute the exact stationary queue length distribution of an M/Ph/m queue, even for its first moment. However, several approximate solutions have been proposed in the literature. Several of these approximations presents an easily implementable closedform formula to evaluate the mean queueing time for a job2 . Possibly inspired by the distributional dependence factor in the Pollaczek-Khintchine formula for the M/GI/1 queue3 , most existing studies seem to concentrate on the influence of the coefficient of variation of the service time distribution. However, a few papers [START_REF] Wolff | The Effect of Service Time Regularity on System Performance[END_REF], [START_REF] Whitt | The effect of variability in the GI/G/s queue[END_REF], and more recently [START_REF] Gupta | The effect of higher moments of job size distribution on the performance of an M/G/s queueing system[END_REF], [START_REF] Brandwajn | Considerations in workload characterization for Parallel Access Volumes[END_REF], report that even the common performance parameters like mean queueing delay depend, sometimes to a large extent, on more than the first two moments of the service time distribution. Hence, it would seem unrealistic to expect an approximation based on only two moments to perform well over a large spectrum of service time distribution. Most of the existing approximations base their solutions on the interpolation and extrapolation of two special cases of the M/Ph/m queue, namely the M/D/m and the M/M/m queues, for which the mean queue length can be easily computed (exactly or approximately). We now present the detailed formula for six of these approximations, often found in textbooks on performance evaluation, to compute the mean queueing time W M/Ph/m . 1) Martin's approximation [START_REF] Martin | System Analysis for Data Transmission[END_REF].

W M/Ph/m P m /µ 1 -ρ . 1 + c 2 B 2m
where P m is the waiting probability for an arriving job estimated by the corresponding probability in an M/M/m queue. 2) Cosmetatos's approximation [START_REF] Cosmetatos | Some Approximate Equilibrium Results for the Multiserver Queue (M/G/r)[END_REF] (also proposed by Björklund and Elldin [START_REF] Björklund | A practical method of calculation for certain types of complex common control systems[END_REF]).

W M/Ph/m c 2 B W M/M/m + (1 -c 2 B )W M/D/m
3) Boxma, Cohen and Huffels referred to as BCH's approximation [START_REF] Boxma | Approximations of the Mean Waiting Time in an M/G/s-Queueing System[END_REF].

W M/Ph/m 1 + c 2 B 2 2W M/D/m W M/M/m 2aW M/D/m + (1 -a)W M/M/m where a = 1 if m = 1 1 m-1 ( c 2 B +1 γ1 -m + 1) if m > 1 and γ 1 = 1 -c 2 B m + 1 + c 2 B m 4
) Tijms's approximation [START_REF] Tijms | Stochastic Modelling and Analysis: A Computational Approach[END_REF].

W M/Ph/m (1 -ρ)γ 1 m + ρ 2 (c 2 B + 1) W M/M/m
where γ 1 is defined as in the BCH approximation. 5) Lee's approximation [START_REF] Lee | Queueing process associated with airline passenger check-in[END_REF].

W M/Ph/m 1 + c 2 B 2 W M/M/m 6) Kimura's approximation [10] W M/Ph/m 1 + c 2 B 2c 2 B W M/M/m + 1-c 2 B W M/D/m
Following guidelines of [START_REF] Bolch | Queueing Networks and Markov Chains[END_REF], we used in the above approximation

W M/D/m 1 2 .W M/M/m 1 + (1 -ρ)(m -1) √ 4 + 5m -2 16ρm .
We implemented other options for W M/D/m (e.g. see eq. (4.17) and (4.29) in [START_REF] Kimura | Approximations for multi-server queues: system interpolations[END_REF]) but they produce virtually identical results in our study. For more details on the rationale behind Martin, Cosmetatos, BCH and Tijms approximations, please refer to [START_REF] Bolch | Queueing Networks and Markov Chains[END_REF], and to [START_REF] Kimura | Approximations for multi-server queues: system interpolations[END_REF] for the others, though their respective domain of applicability is often left unclear.

Overall, little seems to be known about the actual accuracy of these approximations. We briefly review the related state of the art. Not surprisingly, these approximations tend to be fair when the server utilization is low (say less than 0.5). The reason here is that under such condition there is little queue built up in the system, and its performance are mainly driven by the first moment of the service time distribution. Conversely, if the server utilization is high (close to 1), then an M/Ph/m queue is not an adequate model anymore. In such a congestion regime, one needs a queueing model with a finite buffer to assess the system behavior. Existing literature suggests that these approximations will give excellent results when the variability of the service time distribution is low, which means its coefficient of variation ranges from 0 to not much higher than 1. This is quite expected since, by design, most of these approximations rely on the specific solutions of M/D/m and M/M/m queues. Furthermore, statistical distributions with low variability tend to be alike, and thus their expected queueing times tend to be close to those of M/D/m and M/M/m queues. On the other hand, with increasing values of the coefficient of variation of the service time, it is generally known that the approximations accuracy will gradually decrease, although it is generally not stated at what rate. In [START_REF] Bolch | Queueing Networks and Markov Chains[END_REF], the authors compare several approximations introduced above, namely Martin, Cosmetatos, BCH and Tijms, with the results obtained from discrete-events simulation. They consider an M/Ph/m queue with 5 servers, a server utilization rate of 0.7 and values of c B ranging from 0 to 3.2. They conclude that "all approximations are good for c B < 1.4, and that for higher values of c B the approximation due to Cosmetatos is very good and the others are fair."

In this paper, we aim to better characterize the domain of applicability for each of these six approximations. First, we show that their accuracy can be very poor, including for supposedly "easy" examples with a low coefficient of variation of the service time. Second, our results tend to suggest that one can get significantly improved approximate results for the M/Ph/m queue by picking the right approximation according to the actual value of the third moment of the service time, or equivalently, its skewness.

III. NUMERICAL RESULTS, DISCUSSION AND CONCLUSION

A. Experimental protocol

To assess the accuracy of the approximate solutions, and more specifically the influence of the M/Ph/m queue parameters, we proceed as follows. We consider values for the number of servers m ranging from 2 to 12, and for the server utilization ρ ranging from 0.1 to 0.95. For sake of simplicity, we keep the mean service time constant at 1/µ = 1. Taking into account the state-of-the-art, we focus our study on values of the coefficient of variation of the service time c B significantly larger than 1, i.e., ranging from 2 to 10. Finally, unlike most of the existing literature, we pick values for the skewness parameter of the service time s B ranging from c B to 100. Recall that the skewness of a random variable X is defined as follows: s

E[( X-1/µ σ ) 3
] where 1/µ and σ are the mean and the standard deviation, respectively. s is thus the normalized third-order central moment, and it is a measure of the asymmetry of the distribution. The greater the skewness, the longer the tail. Note that a skewness of 0 occurs if the distribution is nearly symmetric. Although in theory negative values are possible for s B , in the case of non-negative distributions and c B > 1 only positive values of s B can occur, meaning that the tail of the distribution is always on the right (see Appendix A). We rely on an easy computable algorithm [START_REF] Bobbio | Matching three moments with minimal acyclic phase type distributions[END_REF] to find a phase type distribution that matches specific values of c B and s B . Then, we compute the exact value of the mean number of jobs in the corresponding M/Ph/m queues using an iterative solution, and evaluate the corresponding relative error for each of the six approximate solutions.

We explore the obtained results along 4 parameters: m, ρ, c B and s B splitting our analysis in two. First we let m and ρ vary while c B and s B are kept constant, and then we do it reverse.

B. Influence of the number of servers and the server utilization

To assess the influence of m and ρ on the approximate solutions accuracy, we consider two distinct examples which we refer to as examples (A) and (B). Their description is given in Table II.

In example (A), the coefficient of variation of the service time is set to 8 while the skewness is equal to 30. Figure 2 shows the degree of accuracy attained by each approximate on this example. It turns out that, with the exception Kimura's solution, virtually all approximations tend to overestimate the queue length, whatever the values of ρ and m. Deviations from the exact solution tend to be higher when ρ is close to 0.5 for small values of m, whereas they peak around ρ equal to 0.7 for larger values of m. Based on this example, the BCH solution with its relative error often less than 10% and at most of 35% significantly outperforms the others whose errors attain and exceed 50%.

Example (B) deals with a smaller value of c B but a larger s B . Figure 3 presents the corresponding results. We observe that, again, all but Kimura's solution tend to overestimate the mean number of jobs in the M/Ph/m queue. Similarly to example (A), it shows that approximate solutions are the most likely to fail when ρ is close to 0.5 for m ranging from 2 to 4. For larger values of m the region of lowest accuracy gradually shifts to higher values of ρ around 0.8. However, unlike example (A), the solution of Kimura yields the best results with a relative error not exceeding 10%.

Based on examples (A) and (B), it seems that when the number of servers is either much larger than the server utilization (say m verifies m > 33ρ -4.66), or conversely much smaller than the server utilization (say m < 50ρ -33), all approximate solutions lead to excellent results. On the other hand, the region where approximate solutions tend to fail occurs for moderate server utilization. When m is low (say less than 4) this region tends to cover server utilizations close to 0.5. As m grows, this region gradually shifts towards server utilizations around 0.8. It is also worth noting that, in this example, BCH's solution always outperforms Martin, Cosmetatos, Tijms and Lee's solution. Finally, it appears that depending on the actual service time distribution (here its c B and s B ), BCH or Kimura's solution emerges as the best choice. In the next section, we study in more detail the influence of the service time distribution.

C. Influence of the coefficient of variation and the skewness

In order to investigate the influence of c B and s B , we now hold constant the number of servers m and the server utilization ρ.

We consider example (C) where m is 4 and ρ is set to 0.7. The corresponding results are illustrated in Figure 4. All solutions except Kimura's approximation yield fair results even for large values of the coefficient of variation as long as the skewness is low (say less than 20). As the skewness increases and reaches 40 or 50, then these solutions lead to highly inaccurate results, unless c B is less than 3. On the other hand, Kimura's solution tends to work best when the other solutions are failing. Except when c B is high and s B is kept low, its deviations from the exact solution are below 25%.

In example (D) we increase the number of servers to m = 8 and we set the server utilization ρ to 0.6. Figure 5 shows the corresponding results. Overall, all the approximate solutions perform better here than in the previous example. Kimura's solution yields excellent results unless c B is high and s B is very low. Conversely, in such a case, all other solutions produce fair results.

From examples of (C) and (D), we observe that if both c B and s B are low (say less than 4 and 30, respectively), any of the six approximate solutions provides accurate results. If the skewness is large (say 10 times greater than the coefficient of variation), but the coefficient of variation is low, our results favor the use of Kimura's solution. Conversely, when the skewness is low, and the coefficient of variation is high, then BCH solution becomes the best option. Finally, when both c B and s B are high, none of the approximations seem to be acceptable.

To conclude, based on our four examples, we recommend the following guidelines: • if m is much larger than ρ (say m > 33ρ -4.66), or conversely much smaller (say m < 50ρ -33), then any approximate solutions yields excellent results; • otherwise if both c B and s B are low (say less than 4 and 30, respectively), any of the six approximate solutions provides accurate results; if s B is large (say 10 times greater than the coefficient of variation) and c B is low, we recommend the use of the Kimura solution; if s B is low and c B is large, we recommend the use of the BCH solution; if both c B and s B are high, none of the approximations seem to be acceptable. The use of these guidelines would ensure that, in almost all our numerous experiments, the discrepancy between the exact and the approximate mean number of jobs in the queue stays within a range of 25% relative error (except obviously when both c B and s B are high and ρ is moderate).

IV. CONCLUSION

The modeling and performance evaluation of systems that process big data is often hampered by difficulties arising due to the high variability of system parameters. In this paper, we assess the potential inaccuracy of several well-established approximations for evaluating the performance of a parallelized device. We focus our study on the case of highly variable and skewed workloads as it is likely to be the case in the context of big data. We use an exact numerical solution to the M/Ph/m queue and an extensive set of experiments (thousands of cases) to evaluate the difference between approximation results and the exact value of the mean queue length. Unlike existing studies, we explore the (in)accuracy of these approximations not only along the commonly explored coefficients of variation but also through the actual value of the third moment of the service time, or equivalently, its skewness.

Our work helps to better characterize the domain of applicability for each of these approximations. First, we show that their accuracy can be very poor, including for supposedly "easy" examples with a low coefficient of variation of the service time. Second, our examples emphasize the important influence of the skewness of the service time distribution on the approximations accuracy whereas none of them account for it. More precisely, our study shows that for 5 out of the 6 tested approximations the results are quite close to each other results (although some differences do exist) and they tend to overestimate the congestion level of the queue. Conversely, the sixth approximate solution is most likely to underestimate the mean length of the queue. Thus, we provide recommendations regarding the choice of the approximation based on the actual set of values for the queue parameters. Using these guidelines would ensure that, in most cases, the relative error of the approximation would remain below 25% relative error.

However, it is important to stress that our results only pertain to a subset of phase type distributions since they were all obtained using a matching algorithm for the first 3 moments of the service time distributions that minimizes the number of phases for the resulting distribution. Thus, a meaningful extension to this work would consist in studying the accuracy of the approximations for phase type distributions that are no more subject to this constraint. Preliminary results (not shown in this paper) seem to indicate that our recommendations still apply in the general case. Another option would be to extend the results of this study to the case of distributions with negative skewness (and hence coefficient of variation less than 1).

Fig. 1 .

 1 Fig. 1. Illustration of an M/Ph/m queue.

Fig. 2 .

 2 Fig. 2. Relative errors of the approximate solutions on the mean number of jobs in the M/Ph/m queue of example (A).

Fig. 3 .

 3 Fig. 3. Relative errors of the approximate solutions on the mean number of jobs in the M/Ph/m queue of example (B).

Fig. 4 .

 4 Fig. 4. Relative errors of the approximate solution on the mean number of jobs in the M/Ph/m queue of example (C).

Fig. 5 .

 5 Fig. 5. Relative errors of the approximate solution on the mean number of jobs in the M/Ph/m queue of example (D).

TABLE II THE

 II M/Ph/m QUEUE PARAMETERS.

	server utilization
	number of servers

A phase type distribution can approximate arbitrarily closely any positive valued distribution[START_REF] Johnson | The denseness of phase distributions[END_REF].

By Little's law[START_REF] Allen | Probability, Statistics, and Queueing Theory with Computer Science Applications[END_REF], the mean queueing time W M/Ph/m and the mean queue length L M/Ph/m are related as follows: L M/Ph/m = λ.W M/Ph/m .

The well-known Pollaczek-Khintchine formula states that the average queueing time in an M/GI/1 queue depends only on the first two moments of the service time.

APPENDIX

A. Relation between coefficient of variation and skewness

Let X be a non-negative random variable. We denote by m i its i-th non-central moment, i.e., m i E(X i ), and by n k its kth normalized moments. By definition, we have:

It is known (e.g., [START_REF] Osogami | Closed form solutions for mapping general distributions to quasi-minimal PH distributions[END_REF]) that for any non-negative distribution we have:

(

Using the definitions of n 2 and n 3 , (1) can be rewritten as:

where s B and c B represent the skewness and the coefficient of variation, respectively. Recall that

and that

. From (2), it follows:

since m 1 and c B are both positive. This simple derivation explains (for a non-negative distribution) why the skewness values can only be positive if we limit the coefficient of variation to values larger than 1.


Stanley Burris

Karen Yeats

Sufficient Conditions for Labelled 0-1 Laws

Keywords: ratio test, labelled structure, zero-one law

If F(x) = e G(x) , where F(x) = P f (n)x n and G(x

This gives an answer to Compton's request in Question 8. 3 [Compton 1987] for an "easily verifiable sufficient condition" to show that an adequate class of structures has a labelled first-order 0-1 law, namely it suffices to show that the labelled component count function is O `nθn ´for some θ ∈ (0, 1). It also provides the means to recursively construct an adequate class of structures with a labelled 0-1 law but not an unlabelled 0-1 law, answering Compton's Question 8.4.

Introduction

Exponentiating a power series can have the effect of smoothing out the behavior of the coefficients. In this paper we look at conditions on the growth of the coefficients of G(x) = g(n)x n , where g(n) ≥ 0, which ensure that f (n -1)/f (n) → ∞, where F(x) = e G(x) . One application of this result is to 0-1 laws, where we find, see Theorem 7, that if the labelled component count function for an adequate class of structures is O(n θn ) for some θ ∈ (0, 1) then the class has a labelled monadic second-order 0-1 law.

Useful notation will be f (n) ≺ g(n) for f (n) eventually less than g(n) and f (n) ∈ RT ∞ for f (n -1)/f (n) → ∞; the notation RT stands for the ratio test. where r n is the unique positive solution to

The Coefficients of e poly

f (n) ∼ F(r n ) r n n • 2πB(r n ) (1
x • G ′ (x) = n,
and B(x

) := x 2 G ′′ (x) + xG ′ (x).
Proof: Theorem X of Hayman [START_REF] Hayman | A generalisation of Stirling's formula[END_REF] shows that F(x) is Hayman-admissible. Then the rest of the claim is an immediate consequence of Corollary II of [START_REF] Hayman | A generalisation of Stirling's formula[END_REF] where the saddle-point method is applied to find the asymptotics of the coefficients of an admissible function. ✷ Corollary 2 For F(x), G(x) as in the above proposition,

(a) f (n) ∈ RT ∞ , (b) f (n) = exp - n log n d 1 + o(1)
.

Proof: Item a follows immediately from Corollary IV of Hayman [START_REF] Hayman | A generalisation of Stirling's formula[END_REF].

For item b one uses r n G ′ (r n ) = n to obtain:

n cdg(d) 1/d r n ≤ n dg(d) 1/d for c > 1 r n = 1 + o(1) n dg(d) 1/d r n n = 1 + o(1) n n dg(d) n/d B(r n ) = 1 + o(1) d 2 g(d) n dg(d) = 1 + o(1) dn G(r n ) = 1 + o(1) g(d)r n d = 1 + o(1) n d F(r n ) = exp n d 1 + o(1) .
Apply these results to [START_REF] Bateman | Monotonicity of partition functions[END_REF]. ✷

Some Technical Lemmas

Now we drop the assumption that G(x) is a polynomial, but keep the requirement

gcd n : g(n) > 0 = 1. (2)
This implies that f (n) ≻ 0. Choose a positive integer L ≥ 2 sufficiently large so

n > L ⇒ [x n ] exp g(1)x + • • • + g(L)x L > 0. (3)
Given ℓ > L with g(ℓ) > 0 let

G 0 (x) := n≥1 g 0 (n)x n := 1≤n≤ℓ g(n)x n F 0 (x) := n≥0 f 0 (n)x n := exp(G 0 (x)) G 1 (x) := n≥1 g 1 (n)x n := n≥ℓ+1 g(n)x n F 1 (x) := n≥0 f 1 (n)x n := exp(G 1 (x)). (4)
Lemma 3 Suppose r ≥ -1 is such that

ng(n) = O f 0 (n + r) . (5)
Then

nf 1 (n) = O f (n + r) .
Proof: In view of (3) and ( 5) we can choose C r such that

ng(n) ≤ C r f 0 (n + r) for n + r ≥ L + 1. (6)
Differentiating (4) gives

nf 1 (n) = n j=ℓ+1 jg(j) • f 1 (n -j) ≤ C r n j=ℓ+1 f 0 (j + r) • f 1 (n -j) by (6) ≤ C r n+r j=0 f 0 (j) • f 1 (n + r -j) = C r f (n + r),
the last line following from F(x) = F 0 (x) • F 1 (x). ✷ Lemma 4 Suppose for every integer r ≥ -1

ng(n) = O f 0 (n + r) . Then f (n -1)/f (n) → ∞. Proof: Since f 0 (n) ∈ RT ∞ by Corollary 2 there is a monotone decreasing function ε(n) such that for any sufficiently large M we have ε(n) > f 0 (n)/f 0 (n -1) for n ≥ M , and ε(n) → 0 as n → ∞.
Thus

f (n) = 0≤j≤n f 0 (j)f 1 (n -j) = 0≤j≤M -1 f 0 (j)f 1 (n -j) + M ≤j≤n f 0 (j)f 1 (n -j) ≤ o f (n -1) + ε(M ) M ≤j≤n f 0 (j -1)f 1 (n -j)
by Lemma 3 and the choice of ε

≤ o f (n -1) + ε(M )f (n -1). Thus lim sup n→∞ f (n) f (n -1) ≤ ε(M ),
and as M can be arbitrarily large it follows that

lim n→∞ f (n) f (n -1) = 0.

✷ 4 Main Result

We are now in a position to prove the main result, making use of

n! = exp n log n • (1 + o(1)) ,
which follows from Stirling's result.

Theorem 5 Suppose F(x) = exp(G(x)) with F(x) = n≥0 f (n)x n , G(x) = n≥1 g(n)
x n , and f (n), g(n) ≥ 0. Suppose also that gcd n : g(n) > 0 = 1 and that for some θ ∈ (0, 1)

g(n) = O n θn n! . Then f (n) ∈ RT ∞ .
Proof: From Corollary 2, for any integer r ≥ -1 and any θ ∈ (0, 1), by choosing ℓ > L such that 1/ℓ < 1 -θ, we have

f 0 (n + r) = exp - (n + r) log(n + r) ℓ 1 + o(1) = exp - n log n ℓ 1 + o(1) ≻ n θn (n -1)! .
Thus ng(n) = O f 0 (n + r) . The Theorem then follows from Lemma 4. ✷

Best Possible Result

The main result is in a natural sense the best possible.

Proposition 6 Suppose t(n) ≥ 0 with gcd n : t(n) > 0 = 1 is such that for every θ ∈ (0, 1)

t(n) = O n θn /n! .
Then there is a sequence g(n) ≥ 0 with gcd n : g(n) > 0 = 1 and g(n

) ≤ t(n) but f (n) / ∈ RT ∞ , where one has F(x) = exp(G(x)). Proof: For θ ∈ (0, 1) let S(θ) = n ≥ 1 : t(n) > n θn /n! .
Then S(θ) is an infinite set.

Let M be such that gcd n ≤ M : t(n) > 0 = 1, and let

g 1 (n) := t(n) if n ≤ M 0 if n > M G 1 (x) := g 1 (n)x n d 1 := deg(G 1 (x)) F 1 (x) := e G1(x) .
For m ≥ 2 we give a recursive procedure to define polynomials G m (x); then letting

d m := deg(G m (x)) F m (x) := e Gm(x) ,
by Proposition 1

f m (n) = exp - n log n d m 1 + o(1) . To define G m+1 (x), having defined G m (x), let h m (n) := 1 n! • n (1-1/2dm)n . Then h m (n) f m (n -1) → ∞ as n → ∞.
Thus we can choose an integer d m+1 > d m such that

d m+1 ∈ S 1 - 1 2d m h m (d m+1 ) > f m (d m+1 -1). This ensures that h m (d m+1 ) ≤ t(d m+1 ). Let G m+1 := G m (x) + h m (d m+1 )x dm+1 . Then f m+1 (d m+1 ) f m+1 (d m+1 -1) ≥ h m (d m+1 ) f m (d m+1 -1) > 1.
Now let G(x) be the nonnegative power series defined by the sequence of polynomials G m (x); and let

F(x) = e G(x) . Then g(n) ≤ t(n) but f (n) / ∈ RT ∞ as f (d m+1 ) f (d m+1 -1) = f m+1 (d m+1 ) f m+1 (d m+1 -1) > 1.
✷ 6 Application to 0-1 laws

A class K of finite relational structures is adequate if it is closed under disjoint union and the extraction of components. One can view the structures as being unlabelled with the component count function p U (n) and the total count function a U (n), both counting up to isomorphism. The corresponding ordinary generating series are

P U (x) := n≥1 p U (n)x n , A U (x) := n≥0 a U (n)x n
connected by the fundamental equation

A U (x) = j≥1 1 -x j -p U (j) . (7)
One can also view the structures as being labelled (in all possible ways) with the count functions p L (n) for the connected members of K, and a L (n) for all members of K. The corresponding exponential generating series are

P L (x) := n≥1 p L (n)x n /n!, A L (x) := n≥0 a L (n)x n /n!
connected by the fundamental equation

A L (x) = e P L (x) . (8)
All references to Compton in this section are to the two papers [START_REF] Compton | A logical approach to asymptotic combinatorics. I. First order properties[END_REF] and [START_REF] Compton | A logical approach to asymptotic combinatorics. II. Monadic second-order properties[END_REF].

Unlabelled 0-1 Laws for Adequate Classes

Let K be an adequate class with unlabelled count functions and ordinary generating functions as described above. Compton showed that if the radius of convergence ρ U of A U (x) is positive then K has an unlabelled 0-1 law (i) iff a U (n) ∈ RT 1 , that is,

a U (n -1) a U (n) → 1 as n → ∞.
K is finitely generated if r = p U (n) < ∞, that is, there are only finitely many connected structures in K. In the finitely generated case the asymptotics for the coefficients a U (n) have long been known to have the simple polynomial form (ii) a U (n) ∼ Cn r-1 (9) provided gcd n : p U (n) > 0 = 1. Item (9) leads to the fact that a U (n) ∈ RT 1 , and hence to an unlabelled 0-1 law. In addition to using this result, Compton notes that the work of Bateman and Erdös [START_REF] Bateman | Monotonicity of partition functions[END_REF] shows that if p U (n) ∈ {0, 1}, for all n, then one has a U (n) ∈ RT 1 .

Both of these results were subsumed in the powerful result of Bell [START_REF] Bell | Sufficient conditions for zero-one laws[END_REF] which says that if p U (n) is polynomially bounded, that is, there is a c such that p U (n) = O(n c ), then a U (n) ∈ RT 1 .

Labelled 0-1 Laws

Compton shows that if ρ L , the radius of convergence of

A L (x), is positive, then K has a labelled 0-1 law iff a L (n -k)/(n -k)! a L (n)/n! → ∞ whenever p L (k) > 0. (10)
In particular it suffices to show that a L (n)/n! ∈ RT ∞ . Compton's method to show that a given adequate class of finite relational structures K has a labelled 0-1 law is to show that its exponential generating function A L (x) = a L (n)x n /n! is Hayman-admissible with an infinite radius of convergence. This guarantees that a L (n)/n! ∈ RT ∞ ([5], Corollary IV). However, as Compton notes, showing that A L (x) is Hayman-admissible can be quite a challenge.

Question 8.3 of [START_REF] Compton | A logical approach to asymptotic combinatorics. I. First order properties[END_REF] first asks if, in the unlabelled case, the result of Bateman and Erdös, namely p U (n) ∈ {0, 1} implies a U (n) ∈ RT 1 , can be extended to the much more general statement that p U (n) = O(n k ) implies a U (n) ∈ RT 1 , yielding an unlabelled 0-1 law. As mentioned earlier, this was proved to be true by Bell. The second part of Question 8.3 asks if there is a simple sufficient condition along similar lines for the labelled case. We can now answer this in the affirmative with a result that is an excellent parallel to Bell's result for unlabelled structures.

Theorem 7 If K is an adequate class of structures with

p L (n) = O n θn
for some θ ∈ (0, 1)

then a L (n)/n! ∈ RT ∞ ,
and consequently K has a labelled monadic second-order 0-1 law.

(i) Given a logic L, K has an unlabelled L 0-1 law means that for any L sentence ϕ, the probability that ϕ holds in K will be either 0 or 1. In [START_REF] Compton | A logical approach to asymptotic combinatorics. I. First order properties[END_REF] Compton worked with first-order logic, in [START_REF] Compton | A logical approach to asymptotic combinatorics. II. Monadic second-order properties[END_REF] with monadic second-order logic. In both papers he simply used the phrases "unlabeled 0-1 law" and "labeled 0-1 law". (ii) This result is usually known as Schur's Theorem [6, 3.15.2]. One can easily find the asymptotics (9) using a partial fraction decomposition of the right side of (7). The labelled case with finitely many components is more difficult-we needed to invoke Hayman's treatise [START_REF] Hayman | A generalisation of Stirling's formula[END_REF] just to obtain the asymptotics for log a L (n)/n! (see Corollary 2). (e) 7.17 Partitions with a Selection Subset p L (n) = 2 n -1.

We can now augment this list by, in each case, coloring the members of K by a fixed set of r colors in all possible ways. This will increase the original p L (n) by a factor of at most r n . This will still give p L (n) = O n n/2 . Furthermore, in each of these colored cases let P be any subset of the connected members, and let K be the closure of P under disjoint union. Each such K has a labelled 0-1 law. Another application of Theorem 7 is to answer Question 4 of [START_REF] Compton | A logical approach to asymptotic combinatorics. I. First order properties[END_REF] by exhibiting an adequate class K such that p L (n) = O n 3n/4 , hence there is a labelled 0-1 law for K; but also such that ρ U ∈ (0, 1), so K does not have an unlabelled 0-1 law.

Let the components of K be the one-element tree T 1 along with rooted trees T 3n of size 3n and height n consisting of a chain C n of n nodes, with an antichain L 2n of 2n nodes (the leaves of the tree) below the least member of the chain; and the chain C n is two-colored while the remaining nodes are uncolored. One can visualize these as brooms with 2-colored handles, see Figure 1.

The number of unlabelled components is given by p U (1) = 1, p U (3n) = 2 n . Thus the radius of convergence of the ordinary generating function of K is ρ U = 3 √ 2. Since this is positive and not 1 it follows from Theorem 5.9(ii) of [START_REF] Compton | A logical approach to asymptotic combinatorics. I. First order properties[END_REF] that K does not have an unlabelled 0-1 law.

For the number p L (3n) of labelled components of size 3n: = (3n) (2/3)(3n) 1+o (1) = O (3n) (3/4)(3n) .

Thus p L (n) = O n 3n/4 , so a L (n)/n! ∈ RT ∞ by Theorem 7, showing that K has a labelled 0-1 law.

Proposition 1 1 F

 11 Given G(x) := g(1)x + • • • + g(d)x d , g(i) ≥ 0, g(d) > 0, with gcd j ≤ d : g(j) > 0 = (x) := n≥0 f (n)x n = e G(x) ,the function F(x) is Hayman-admissible. Thus

  ) 1365-8050 c 2008 Discrete Mathematics and Theoretical Computer Science (DMTCS), Nancy, France

Fig. 1 :

 1 Fig. 1: Brooms with two-colored handles

p≤ 2 n

 2 L (3n) ≤ 2 n 3n n n! (3n) n exp n log n • (1 + o(1)) = exp 2n log n • (1 + o(1))