John W Vasquez
email: jwvasque@laas.fr

Louise Travé-Massuyès
email: louise@laas.fr

Audine Subias
email: subias@laas.fr

Fernando Jimenez
email: fjimenez@uniandes.edu.co

Carlos Agudelo
email: carlos.agudelo@ecopetrol.com.co

Chronicle based alarm management in startup and shutdown stages

The transitions between operational modes (startup/shutdown) in chemical processes generate alarm floods and cause critical alarm saturation. We propose in this paper an approach of alarm management based on a diagnosis process. This diagnosis step relies on situation recognition to provide to the operators relevant information on the failures inducing the alarms flows. The situation recognition is based on chronicle recognition. We propose to use the information issued from the modeling of the system to generate temporal runs from which the chronicles are extracted. An illustrative example in the field of petrochemical plants ends the article.

Introduction

The petrochemical industries losses have been estimated at 20 billion dollars only in the U.S. each year, and the AEM (Abnormal Events Management) has been classified as a problem that needs to be solved. Hence the alarm management is one of the aspects of great interest in the safety planning for the different plants. In the process state transitions such as startup and shutdown stages the alarm flood increases and it generates critical conditions in which the operator does not respond efficiently, then a dynamic alarm management is required [START_REF] Zhao | A dynamic alarm management strategy for chemical process transitions[END_REF]. Currently, many fault detection and diagnosis techniques for multimode processes have been proposed; however, these techniques cannot indicate fundamental faults in the basic alarm system [START_REF] Thornhill | A combined analysis of plant connectivity and alarm logs to reduce the number of alerts in an automation system[END_REF], in the other hand the technical report "Advance Alarm System Requirements" EPRI (The Electric Power Research Institute) suggests a cause-consequence and event-based processing. In this perspective, diagnosis approaches based on complex events processing or situation recognition are interesting issues. Therefore, in this paper, a dynamic alarm management strategy is proposed in order to deal with alarm floods happening during transitions of chemical processes. This approach relies on situations recognition (i.e. chronicle recognition). As, the efficiency of alarm management approaches depends on the operator expertise and process knowledge, our final objective is to develop a diagnosis approach as a decision tool for operators. The paper is divided into 6 sections. Section 2 gives an overview on the relevant literature. The section 3 concerns the modeling of the system. The section 4 is about the chronicle principle and the temporal runs used for the chronicle design. The section 5 is devoted to the chronicle generation. Finally , an illustrative application on real data from a petrochemical plant is given section 6.

State of the art: Alarm management

Alarm management has recently focused the attention of many researchers in themes such as:

Alarm historian visualization and analysis: A combined analysis of plant connectivity and alarm logs to reduce the number of alerts in an automation system was presented by [START_REF] Izadi | Graphical tools for routine assessment of industrial alarm systems[END_REF]; the aim of the work presented is to reduce the number of alerts presented to the operator. If alarms are related to one another, those alarms should be grouped and presented as one alarm problem. Graphical tools for routine assessment of industrial alarm systems was proposed by [START_REF] Takai | Use of event correlation analysis to reduce number of alarms[END_REF], they presented two new alarm data visualization tools for the performance evaluation of the alarm systems, known as the high density alarm plot (HDAP) and the alarm similarity color map (ASCM). Event correlation analysis and two-layer cause-effect model were used to reduce the number of alarms in [START_REF] Ge | Multimode process monitoring based on bayesian method[END_REF]. A Bayesian method has been introduced for multimode process monitoring in [START_REF] Noda | Evaluation of plant alarm systems by behavior simulation using a virtual subject[END_REF]. This type of techniques helps us to recognize alarm chattering, grouping many alarms or estimate the alarm limits in transition stages, but the time and the procedure actions are not included.

Process data-based alarm system analysis and rationalization: The evaluation of plant alarm systems by behavior simulation using a virtual subject was proposed by [START_REF] Yang | Improved correlation analysis and visualization of industrial alarm data[END_REF]. [START_REF] Shook | A framework for optimal design of alarm systems[END_REF] introduced a technique for optimal design of alarm limits by analyzing the correlation between process variables and alarm variables. In 2009 a framework based on the receiver operating characteristic (ROC) curve was proposed to optimally design alarm limits, filters, dead bands, and delay timers; this work was presented in [START_REF] Oktem | Dynamic risk analysis using alarm databases to improve process safety and product quality: Part ii bayesian analysis[END_REF] and a dynamic risk analysis methodology that uses alarm databases to improve process safety and product quality was presented in [START_REF] Liu | Non stationary fault detection and diagnosis for multimode processes[END_REF]. In [START_REF] Boang | Fault diagnosis strategy for startup process based on standard operating procedures[END_REF] the Gaussian mixture model was employed to extract a series of operating modes from the historical process data and then the local statistic and its normalized contribution chart were derived for detecting abnormalities early and for isolating faulty variables. We can see that the use of virtual subjects could be applied to probe the alarm system and using historical information about the alarm behavior for detecting abnormalities. The problem is presented when the simulation requires a lot time to probe the totally of scenarios and when we have new plants that do not contain information about historical data.

Plant connectivity and process variable causality analysis (causal methods): In the literature, transition monitoring of chemical processes has been reported by many researchers. In [START_REF] Vedam | A framework for managing transitions in chemical plants[END_REF] was presented a fault diagnosis strategy for startup process based on standard operating procedures, this approach proposes behavior observer combined with dynamic PCA (Principal Component Analysis) to estimate process faults and operator errors at the same time, and in [START_REF] Travé-Massuyès | Learning chronicles signing multiple scenario instances[END_REF] was presented a framework for managing transitions in chemical plants where a trend analysis-based approach for locating and characterizing the modes and transitions in historical data is proposed. Finally, in [?] a hybrid model-based framework was used for alarm anticipation where the user can prepare for the possibility of a single alarm occurrence. For the transition monitoring, these types of techniques are the most used in industrial processes and the hybrid model based framework could be a good representation of our system. We can observe that a causal model allows identify the root of the failures and check the correct evolution in a transitional stage. Our proposal is closer to the third type of approach and seeks to exploit the causal relationships as presented in the next sections.

3 Representation of the system

Hybrid Causal Model

The hybrid system is represented by an extended transition system [?], whose discrete states represent the different modes of operation for which the continuous dynamics are characterized by a qualitative domain. Formally, a hybrid causal system is defined as a tuple:

Γ = (ϑ, D, Conf, T r, E, CSD, Init) (1)
Where

• ϑ = {v i } is a set of continuous process variables which are function of time t. • D is a set of discrete variables. D = Q ∪ K ∪ V Q . Q
is a set of states q i of the transition system which represent the system operation modes. The set of auxiliary discrete variables K = {K i , i = 1, ...n c } represents the system configuration in each mode q i as defined below by Conf (q i ). V Q is a set of qualitative variables whose values are obtained from the behavior of each continuous variable v i .

• Conf (q i ): Q → ⊗ i D(K i )
where ⊗ is the Cartesian product and D(K i ) is the domain of K i ∈ K that provides the configuration associated to the mode. i.e. the modes of the underlying multimode components (typically, a valve has two normal modes, opened and closed) • E = Σ∪Σ c is a finite set of event types noted σ, where:

-Σ is the set of event type associated to the procedure actions in a startup or shutdown stages. -Σ c is the set of event type associated to the behavior of the continuous process variables. • T r : Q× Σ → Q is the transition function. The transition from mode q i to mode q j with associated event σ is noted (q i , σ, q j ) or q i σ -→ q j . We assume that the model is deterministic, without loss of generality i.e. whenever q i σ -→ q j and q i σ -→ q k then q j = q k for each (q i , q j , q k ) ∈ Q 3 and each σ ∈ Σ.

• CSD

⊇ i CSD i is the Causal System Description or the causal model used to represent the constraints underlying in the continuous dynamic of the hybrid system. Every CSD i associated to a mode q i , is given by a graph (G c = ϑ ∪ K, I). I is the set of influences where there is an edge e(v i , v j ) ∈ I from v i ∈ ϑ to v j ∈ ϑ if the variable v i influences variable v j . Then, the vertices represent the variables and the edges represent the influences between variables and for each edge exists an association with a component in the system. The set of components is noted as COM P .

• Init is the initial condition of the hybrid system,

Qualitative abstraction of continuous behavior

In each mode of operation, variables evolve according to the corresponding dynamics. This evolution is represented with qualitative values. The domain

D(V i ) of a qualitative variable V i ∈ V Q is obtained through the function f qual : D(v i ) → D(V i
) that maps the continuous values of variable v i to ranges defined by limit values (High H i and Low L i ).

f (v i ) qual =              V H i if v i ≥ H i ∧ dvi dt > 0 V M i if v i < H i ∧ dvi dt < 0 ∨ v i ≥ L i ∧ dvi dt > 0 V L i if v i < L i ∧ dvi dt < 0 (2)
dvi dt > 0 represents that the continuous variable v i is increasing and dvi dt < 0 that it is decreasing. The behavior of these qualitative variables is represented in Figure 1. by the graph

G Vi = (V Q , Σ c , γ) where V Q is the set of the possible qual- itative states (V L i : Low, V M i : M edium, V H i :
High) of the continuous variable v i , Σ c is the finite set of the events associate to the transitions and γ : V Q × Σ c → V Q is the transition function. The corresponding event generator is

V Q →σ f V Q →σ : V Q × γ(V Q , Σ c ) → Σ c ∀V i ∈ V Q , (V n i , V m i ) →        l + (v i ) if V L i → V M i l -(v i ) if V M i → V L i h + (v i ) if V M i → V H i h -(v i ) if V H i → V M i V n i , V m i ∈ {V L i , V M i , V H i } (3) Σ c = vi∈ϑ {l + (v i ), l -(v i ), h + (v i ), h -(v i )} (4)

Automatic derivation of the causal model

To obtain the causal model of a system in a given operating mode implies to collect the equations that represent the behavior of the system in this mode. The theory of causal ordering issued from the Qualitative Reasoning community can be well applied to obtain automatically the causal structure associated to a set of equations. Now, associating activation conditions to the equations extend the causal ordering to systems with several operating modes [START_REF] Guerraz | Chronicles construction starting from the fault model of the system to diagnose[END_REF]. Then these activation conditions can be related in the influences of the resulting causal graph.The proposed algorithm, implemented in the Causalito software makes use of conditions that avoid recomputing a totally new perfect matching for every operating mode, thus reducing the computational cost. In this work, the Causal System Description is given by CSD = (ϑ, I), where each influence I is labeled with:

• An activation condition indicating the modes in which it is active (or no label if it is active in all modes), • The corresponding equation, • The component whose behavior is expressed by the equation. In the follow section we expose the principle of the chronicle generation where concepts such as event, chronicle and temporal run are described.

Chronicles

Events and chronicles

Let us consider time as a linearly ordered discrete set of instants. The occurrence of different events in time represents the system dynamics and a model can be determined to diagnose the correct evolution. An event is defined as a pair (σ i , t i ), where σ i ∈ E is an event type and t i is a variable of integer type called the event date. We define E as the set of all event types and a temporal sequence on E is an ordered set of events denoted S = (σ i , t i ) j with j ∈ N l where l is the size of the temporal sequence S and N l is a finite set of linearly ordered time points of cardinal l. A chronicle is a triplet C = (ξ, C T , G) such that ξ ⊆ E, C T is the set of temporal constraints. G = (N, It) is a directed graph where N represent event types of E and the arcs It represent the relationship between events σ ∈ E, if the event σ 1 occurs t time units after σ 2 , then it exists a directed link from σ 1 to σ 2 associated with a time constraint. Considering the two events (σ i , t i ) and (σ j , t j ), we define the time interval as the pair τ ij = [t -, t + ], τ ij ∈ C T corresponding to the lower and upper bounds on the temporal distance between the two event dates t i and t j [21]. The idea of our proposal is to design the chronicles from the hybrid causal model of the system. Indeed the evolution of the system can be captured with temporal runs from which chronicles can be learn (See Figure 2). More precisely, the system initiates in the state q 0 and it evolves according to the transitions resulting from the events defined by the procedure actions for specific scenarios (startup/shutdown). For a given system modes q i ∈ Q, the associated CSD i is used to generate the set of event types corresponding to the evolution of the continuous process variables. A run is defined by a sequence of event types α 1 , α 2 , ....α n where α i ∈ E generated for each scenario using the startup/shutdown procedures. These runs with time constraints permit to construct the chronicle database of the system. In this preliminary approach, time constraints are obtained by simulation.

Temporal runs

We denote a temporal run as R, T where R is a run and T is the time graph of the run that includes the time constraints C T between each pair of time points where must occurs the events type. Figure 3 gives time graph examples and the possible composition of time graphs. In our approach the Figure 3: Time graphs example runs are issued from the system evolution from one operation mode to another. The interleaved sequence of event types α 1 , α 2 , . . . α n represents the procedure actions and the behavior evolution of the process variables. The time constraints between each pair of event types are determined by simulation of the continuous behavior for each process variable, responding to the procedure actions. (5) The chronicle database used for diagnosis is composed by the entries of all the matrices {CAr k }. This chronicle database is submitted to a chronicle recognition system that identifies in an observable flow of events all the possible matching with the set of chronicles from which the situation (normal or faulty) can be assessed.

Generation of Chronicles

CAr k = O 1 O 2 O 3 O 4 . . . . . . O j N f 1 f 2 . . . . . . f n          C k 01 C k 11 C k 21 . . . . . . C k i1 C k 02 C k 12 C k 22 . . . . . . C k i2 C k 03 C k 13 C k 23 . . . . . . C k i3 C k 04 C k 14 C

Chronicle learning

As explained previously when the system changes mode of operation, a set of event types occurs forming a run R. As this evolution is due to procedure actions. Not only a unique temporal run can occur. Hence, we need to set up the maximal number of temporal runs that it could occur in each scenario represented in the matrix [START_REF] Takai | Use of event correlation analysis to reduce number of alarms[END_REF]. To obtain the chronicle in each scenario is necessary to obtain the larger time graph with as many event types and with the minimal values of the constraints. [19] proposes to determine the chronicles from the temporal runs. They define a partial order relation between two temporal runs as R, T ≤ R , T when the set of event types in R is a subset of event type in R and the time graphs T and T are related by T T determining the result graph where exists a unique equivalent constraint that is the minimal. The relation expresses that the set of constraints in the time graph T is a subset of constraints in T , C T (t, t ) ⊆ C T (t, t ). Therefore, we apply the composition (see Figure 3) between the time graphs in order to merge the constraints obtaining the larger and constrained time graph that represents the chronicle in that scenario. Figure 4 gives an example of a chronicle generation from a maximal temporal run. In the next section a case study is presented in Figure 4: Chronicle example which the chronicle generation from the temporal runs is illustrated.

6 Case study 6.1 HTG (Hydrostatic Tank Gauging) system

In the Cartagena Refinery currently are being implemented news units and elements. In the startup stage they will need a tool to help the operator to recognize dangerous conditions. We will analyze the startup and shutdown stages in the unit of water injection. This process is a HTG (Hydrostatic Tank Gauging) system composed by the following components: one tank (T K), two normally closed valves (V 1 and V 2), one pump (P u), a level sensor (LT ), a pressure sensor (P T ), inflow sensor (F T 1 ) and an outflow sensor (F T 2 ), see Figure 5. Assuming this system as a hybrid causal model, the underlying discrete event system and the different process operation modes are described in Figure 6 where we can see a possible correct evolution for the startup procedure. The events V 1 c,o , V 2 c,o represents that the valves V 1,V 2 move from the state closed to the state opened, the events V 1 o,c ,V 2 o,c represents on the contrary the valves moving from the state opened to the state closed. The event P u f -n indicates that the pump P u is turned on and the event P u n-f indicates that the pump P u is turned off.

Identification of causal relationships

The level (L) in the tank is related to the weight (m) of the liquid inside, its density (ρ) and the tank area (A). The density (ρ) is the relationship of the pressures (P med ,P inf ) in separated points (h). Based on the global material balance, we define that the input flow is equal to the outlet flow. Then, the variation of the weight (dm(t)/dt) in the tank is proportional to the difference between the inflow (Q iT K ) and the outflow (Q oV 2 ). The differential pressure in the pump and in V 2 are specified as ∆ P P u and ∆ P V 2 . The outlet pressure in the pump (P o) is related with the outlet flow tank (Q oT K ), the revolutions per minute in the pump (RP M P u ), his capacity (C) and the radio of the outlet pipe (r). The outflow (Q oV 2 ) and inflow (Q iT K ) control are related to the percentage aperture of the valves V 1 (LV 1) and V 2 (LV 2) and differential pressures (∆P V 1 ,∆P V 2 ). In Figure 7 we can see the CSD of the system in the modes q 1 , q 5 and q 7 . For example, the mode q 1 activates the influence of Q iT K to L. The mode q 5 activates the influence of Q iT K to L and the influence of L to P o and finally the mode q 7 activates the influence of Q iT K to L, L to P o and P o to QoV 2. 7: CSD in the modes q 1 , q 5 and q 7

Event identification

One of the most important steps for fault diagnosis based on chronicle recognition is to determine the set of events that can carry the system to a failure. Each situation pattern (normal or abnormal) is a set of events and temporal constraints between them; then a situation model may also specify events to be generated and actions to be triggered as a result of the situation occurrence. For a startup procedure in the example process, the set of event types Σ that represent the procedure actions is:

Σ = {V 1 c,o , V 2 c,o , P u f -n , V 1 o,c , V 2 o,c , P u n-f } (6)
According to the causal graphs associated to the modes involved in the sequence of procedure actions (i.e q 1 , q 5 and q 7 indicated by red arrows on Figure 6), the event types of Σ c correspond to the behavior of the variables L,P o and Q oV 2 .

Σ c = {l + (L) , l - (L) , h + (L) , h - (L) , l + (P o) , l - (P o) , h + (P o) , h - (P o) , l + (Q oV 2 ) , l - (Q oV 2 ) , h + (Q oV 2 ) , h - (Q oV 2 ) } (7)
From the startup/shutdown procedures the different temporal runs are determined and these temporal runs are related to the normal and abnormal situations. The chronicle resulting from a normal startup procedure is presented in Figure 8. The model system was developed in Matlab including A preliminary method for alarm management based on automatically learned chronicles has been proposed. The proposal is based on a hybrid causal model of the system and a chronicle based approach for diagnosis. An illustrative example of an hydrostatic tank gauging has been considered to introduce the main concepts of the approach. In this paper the design of the temporal constraints of the chronicles were performed from simulation results, but further research aim to generate the chronicles from the model of the system. Learning approaches are currently considered for acquiring the chronicle base directly from the sequences of events representing the situations. For this propose the algorithm HC-DAM (Heuristic Chronicle Discovery Algorithm Modified [21]) may be used. The use of HIL (Hardware in the loop) to simulate and validate the proposal is also in our prospects.

Acknowledge

The ECOPETROL -ICP engineers Jorge Prada, Francisco Cala and Gladys Valderrama help us to develop and validate the simulations.

Figure 1 :

 1 Figure 1: Qualitative values of the process variables defined by the abstraction function f V Q →σ

Figure 2 :

 2 Figure 2: Principle of chronicle generation

5. 1

 1 Chronicle databaseAn industrial or complex process P r is composed of different areas P r = {Ar 1 , Ar 2 , ...Ar n } where each area Ar k has different operational modes such as startup, shutdown, slow march, fast march, etc. The set CAr k of chronicles C k ij for each area Ar k is presented in the matrix below, where the rows represent the operating modes (i.e. O 1 : Startup, O 2 : Shutdown, O 3 : Startup type , O 4 : Startup type , etc) and the columns the different faults.

  k 24 . . . . . . C k i4 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Figure 5 :

 5 Figure 5: Process diagram

Figure 6 :

 6 Figure 6: Underlying DES of the HGT system

Figure

  Figure7: CSD in the modes q 1 , q 5 and q 7

Figure 8 :

 8 Figure 8: Chronicle C 01 for normal behavior startupthe injection water process area. The continuous behavior is related to the evolution of the level L, outlet pump pressure P o and the outlet flow Q oV 2 in the system. The discrete evolution is related to the event evolution of the procedures in the startup and shutdown stages. From the different failure modes of the process, the dynamic behavior of the variables is shown with a detection for the possible process states, including the normal procedure without failure. The simulation includes 3 types of startup procedures (OK, f ail 1 and f ail 2 ) with 4 types of fault modes (V 1 , V 2 , P ump and Drain open ) and 3 types of Shutdown procedure (OK, N on -actived and F ail). The evolution of the continuous variables in the startup procedure without failure is shown in Figure9. The events are generated by the program through the evolution of the differential equations, the

Figure 9 :

 9 Figure 9: Normal behavior in startup procedure without failure. Blue: Level, Green:Pressure, Red: ouletflow


Tom Obry

Audine Subias

Louise Travé-Massuyès

A Learning Algorithm for Episodes

Sequences of events describing the behavior and actions of agents or systems can be collected in several domains. An episode is a collection of events that occur in a given partial order. By performing a recognition of recurrent episodes in several sequences and comparing them, it is possible to determine a pattern common to all the sequences. In this paper, we propose an approach to recognize episodes that are common in a set of event sequences. The method described is applied to the automotive domain for learning diagnosis procedures.

Introduction

In the automotive domain, as vehicles become increasingly complex, it becomes difficult to efficiently train garage mechanics in front of the ever growing technology. Faults become more and more difficult to diagnose and training courses multiply to manage new failures. The purpose of this work is to improve the daily work of garage mechanics by providing them with assistance in the diagnosis and repair tasks.

A diagnosis session corresponds to a set of tests and actions performed by the mechanics via a diagnostic tool performed on a vehicle with a breakdown. Each diagnosis session performed is saved in the form of a trace. The goal is to compare all of these diagnosis traces corresponding to the same symptoms to extract the common actions. The aim of this work is to determine the most common action plan that successfully diagnoses a breakdown or that carries out a repair. We propose an approach based discovery and learning patterns [START_REF] Agrawal | Fast algorithms for mining association rules[END_REF] [START_REF] Mitsa | Temporal data mining[END_REF]. It is inspired by [3] [4] that proposes an algorithm to learn chronicles. Chronicles are temporal patterns based on a set of events and temporal constraints relating their occurence dates. They are used to represent the dynamic behavior of complex systems in an event-based fashion and have found several applications in the supervision and diagnostic field [5; 3; 6]. In our case, the occurence dates are not recorded in the traces that are given as input sequences for the learning algorithm. The reason is that actions performed by the mechanics are not constrained over time. This is why we turned to a particular type of chronicles called "episodes" [START_REF] Verkalo | Discovery of frequent episodes in event sequences[END_REF] [START_REF] Amphawan | Mining top-k regular episodes from sensor streams[END_REF]. Episodes just specify the precedence order of events. This paper is organized as followed. Section 2 introduces the notion of episode. In the section 3, a new algorithm to learn episodes (LAE) is described. Section 4 presents the test and results. Section 5 presents a standard algorithm for episode learning and a comparison with the new algorithm. Finally, section 6 concludes with the perspectives of this work. 1

Basic concepts

We consider the input as a sequence of events, where each event has an associated date of occurrence. Given a set E of event types, an event is a pair (e i , t i ), where e i ∈ E is an event type and t is an integer, the (occurrence) date of the event.

Definition 1. An event sequence s on E is a triple (s, T s , T e ) s =< (e 1 , t 1 ), (e 2 , t 2 ), ..., (e n , t n ) > (1)
is an ordered sequence of events such that e i ∈ E for all i = 1, ..., n, and t i ≤ t i+1 for all i = 1, ..., n -1. Further on, T s and T e are integers: T s is called the starting time and T e the ending time, and T s ≤ t i ≤ T e for all i = 1, ..., n [START_REF] Verkalo | Discovery of frequent episodes in event sequences[END_REF].

Episodes are defined as a collection of events that occur frequently in a sequence. An episode is a set of partially or totally ordered event types [START_REF] Verkalo | Discovery of frequent episodes in event sequences[END_REF]. Definition 2. An episode ψ = (V, ≤, g) is a set of nodes V , an ordered partial ≤ on V and a mapping g : V → E which associates a node with an event type.

The interpretation of an episode means that the events of g(V ) must appear in the order described by ≤. The difference with chronicles is at the level of the temporal bounds: the intervals in episodes are not bounded. This makes it possible to consider only the notion of precedence between the events:

t j -t i ∈ [0, +∞] or t j -t i ∈ [-∞, 0]. Definition 3. An episode β = (V , ≤ , g ) is a sub-episode of α = (V, ≤, g), denoted β ≤ α if there exists an injective mapping f : V → V such that g (v) = g(f (v)) for all v ∈ V and for all v, w ∈ V with v ≤ w also f (v) ≤ g(w).
An episode α is a super-episode of β if and only if β ≤ α. There are two classes of episodes: serial episodes and parallel episodes. Definition 4. An episode α is serial if the relation ≤ is a total order (i.e., x ≤ y or y ≤ x for all x, y ∈ V ). If the partial order ≤ is trivial (i.e., x y for all x, y ∈ V such that x = y), the episode is parallel [START_REF] Mannila | Discovering generalized episodes using minimal occurrences[END_REF]. Figure 1 shows two examples of episodes. The episode on the left is a serial episode because the order is total between A, B and C. The other episode is parallel. Indeed, there is no constraint on the partial order on the event types A, B or C [START_REF] Verkalo | Discovery of frequent episodes in event sequences[END_REF]. Episodes issued from the diagnostic tool traces are parallel episodes because the order between the event types is not total. This means that the order of the actions for diagnose a vehicle is not total.

Learning Algorithm for Episodes (LAE)

We present in this section, a new algorithm which can compute a set of parallel episodes from a set of input sequences. The solution includes events that are common to all input sequences. We detail the different steps that make up the algorithm. For a better comprehension of the steps, we use 3 sequences, s 1 ∈ S, s 2 ∈ S and s 3 ∈ S (where S is the set of the input sequences): The sequences are stored in a table. Each cell contains 2 other cells, where the first contains the event types and the second, the dates of occurrence.

s 1 =< (A,

Searching for common events

The first step consists in sorting the events that are common and those that are not. With these common events, events that are not in the initial sequences are removed. The k th event e k of the i th sequence s i is compared with the j th event e j of the collection of common events Com[j]. When the two events compared are identical, the algorithm indexes the event in the first cell and the date of occurrence in the second cell on s i . 2 summarizes, for each sequence, the event types that are common in all sequences and their occurrence positions in the sequence. All the events located in s 1 are in Com [START_REF] Agrawal | Fast algorithms for mining association rules[END_REF] because they are present in all the sequences. The event type E is not present in the s 2 and s 3 because E is not present in all the sequences.

Computing the number of occurrences of an event type in n sequences and computing the pairs of events

In this step, the number of occurrences of each event type is found in the n input sequences. These information are stored in a table called Set. It is composed of k columns and i rows, where k corresponds to the event type and i, the input sequence number. Each cell is composed of 2 more cells. The first contains the frequency of e k located in the i th sequence. The second cell contains a vector of "0" and "1". "1" indicates the position of the event e k in the sequence i.

A B C D 2 [100010] 1 [010000] 1 [001000] 1 [000100] 2 [100001] 1 [010000] 1 [001000] 1 [000100] 2 [001010] 1 [100000] 1 [010000] 1 [000100]
Table 3: Synthesis of the information in Set Table 3 shows the event types common in Set, their occurrence positions for each sequence and the position of each event in each sequence. To determine the maximal frequency of each event, noted F req(e k ), the minimal number of occurrences of each event, noted min(Occ(e k )), for each sequence s is required.

F req(e k ) = min(Occ(e k ) ∈ s, s ∈ S, ∀s) (5)
Each frequency is stored in a table F req, under the format (e k , F req(e k )). All the pairs and their maximal frequencies in S are presented in the table 5. In this case, the maximal frequency for pairs (A, B), (A, C) and (A, D) is 2 and for pairs (B, C), (B, D) and (C, D), the value of the maximal frequency is 1.

A B C D 2 1 1 1

Determination of the direction of the pairs

We are going to explain, for each pair of event types, a relation of unique precedence, that is to say present on the n input sequences. In other words, it is a matter of determining the direction of a pair of events. Let's take our 3 examples: For the pair (A, B), A is before and after B in s 1 and s 2 .

s 1 =< (A,
The label "↔" is allocated to the pair (A,B) to indicate that the pair (A, B) can occur in any order, i.e. A can precede B or B can precede A. For example, the pair A ↔ B means it is possible to go to the event B from A and return to the event A. Pairs (B, C), (B, D) and (C, D) keep always the same direction in all sequences. The label is given for those two pairs of events "→" because C always follows B, D is always after B and D always occurs after C. To summarize, the label "→" can only be assigned if A is followed by B in all the sequences (conversely, "←" for B is followed by A in the all the sequences). The label "↔" is assigned to a pair:

• if there is a label "→" and a label "←" in at least one sequence. It indicates that both directions are possible.

• if the two directions are possible in the same sequence.

The table Set is organized by events in columns indexed k.

Every event e k gives rise to a macro-column with two subcolumns indexed k 1 and k 2 . The lines indexed i correspond to the input sequences. The vectors of "0" and "1", located in table Set are studied because they contain the occurrence positions of the events. Set(i, k 1 ) provides the number of occurence of the event of column k in the sequence i and Set(i, k 2 ) is the vector from which we can retrieve the event position in the sequence. From this information, we can obtain the direction of every pair in every sequence. To save the direction of the pairs, the variable Inter is defined and can take several values:

• Inter i (k, l) =+1 if e l follows e k in the i th sequence,
• Inter i (k, l) =-1 if e k follows e l in the i th sequence,

• Inter i (k, l) = 0 if e k precedes and follows e l in the i th sequence.

To find the direction of the pairs in the n sequences, we proceed as follows:

• if n i=1
Inter i (k, l) = n, the label assigned to the pair (e k , e l ) is "→",

• if n i=1 Inter i (k, l) = -n
, the label assigned to the pair (e k , e l ) is "→",

• otherwise, the label assigned to the pair (e k , e l ) is "↔".

Test of consecutivity

We aim to reduce the set of solutions by adding the notion of consecutiveness. If two events are always consecutive on the n sequences (i.e there is no other event between them) then we can say that the pair is consecutive. In this case, a "No_Event" constraint between the two events is set [START_REF]Discovering chronicles with numerical time constraints from alarm logs for monitoring dynamic systems[END_REF]. The vectors of the occurrence positions Set(i, k 2 ) of the first pair events e k and e l are studied. When the first "1" of the e k vector is detected, the value of the following column of the vector of the event e l is studied: if it is a "1", then the pair studied is consecutive at least once in one sequence and the algorithm continues the test with the other sequences and with the others pairs. When all pairs are calculated, they are assembled to be interpreted with a contraints graph. A contraints graph represents the result with a graphic way. Each node represents an event and the index next to the event indicates the maximal frequency of each event. The arc corresponds of the no bounded time contraints. A dashed line represent the specific relation "No_Event" between a pair of events. Arrows on the arcs model the direction of each pair.

Tests and results

We applied the LAE algorithm with the ACTIA's traces. The actual traces were not available during the tests, they were simulated on a diagnostic tool. The diagnostic tool allows the mechanic to carry out his tests and measurements on the vehicle in real time. The mechanic selects actions displayed on the screen. The tool works by changing screens. The users can interact via the tool to test the various functions of the vehicle. The extract from figure 4 represents the beginning of a simulated trace on the diagnostic tool. The first accessible screen offers a choice on the brand of the vehicle to be diagnosed, "Peugeot" or "Citröen". This example was carried out on a vehicle "Citröen" and the model was a "C4 Picasso". The first line corresponds to the selection of the "Citröen" button of the diagnostic tool. A change of screen displays a list of vehicles of the mark to be selected. The second line indicates that the mechanic was on the "Citröen" vehicle selection menu and the third line shows he selected the "C4PB78" button, which corresponds to the "C4 Picasso" model, etc. The simulated sequences have an average of 65 events and there are 4. The result consists of 78 pairs including 13 elements in common. Note that in this case of study, two actions (i.e. events) without preference of order cannot be done sequentially in one or the other order equivalently. The generated result should be filtered to take this feature into account. The algorithm is able to consume traces of ACTIA and generate a satisfactory result from a temporal and qualitative point of view. The execution time is 0.7 seconds with ACTIA traces. The tests were carried out on an ASUS PC, under windows 7 with an Inter Core i7 processor -2670QM, 2.2GHz, 6GB of RAM. The work carried out is considered as a proof of concept. Nevertheless, this work need to be continued by a deeper analysis of the results. Indeed, the constraints graph represents the raw episode of the n input sequences. As the result was not confronted with the mechanics, the pertinence of the episode was not evaluated.

5 A standard algorithm for episode learning

Presentation of the algorithm

In this section, we present an existing algorithm that allows learning episodes in n sequences [START_REF] Verkalo | Discovery of frequent episodes in event sequences[END_REF]. This algorithm allows to extract the different serial episodes and parallel episodes that are located in a sequence s. The sequences provided by the diagnostic tool being parallel, we only present the part dealing with this type of episodes. The algorithm recognizes parallel episodes in a sequence s. This step creates the database that contains all the episodes of the s sequence. The algorithm uses two windows w = (w, t s , t s + win) et w = (w , t s + 1, t s + win + 1) where w and w are the sequences contained mutually in w and w , t s is the starting time of the window and win is the dimension of a window such that width(w) = win [START_REF] Verkalo | Discovery of frequent episodes in event sequences[END_REF]. The sequences contained in w et w will therefore be similar. After the recognition of the episode in w, an update of the data structure is performed to get the window to move to the w episode in the same sequence s. The recognition of a candidate parallel episode α is done with a counter α.event_count. This counter indicates how much events of α are present in the window. When α.event_count is equal to the length of α, the totality of α is present in the window. α.f req_count is incremented by 1. That counter represents the number of windows where α is contained entirely. Candidates episodes are indexed by the number of events of each type they contain. All episodes that contain a type A events will be in the contain(A, a).

An episode α is stored in a table of events sorted alphabetically. Each event of α is indexed between brackets. For example, an episode β with events types B, D, D, E is represented by a table

β with β[1] = B, β[2] = D, β[3] = D, β[4] = E.
Collections of episodes are sorted alphabetically in a table F , where the i th episode of the collection is represented by F [i]. All the episodes that share the first event type are consecutive in the episode collection. A maximum sequence of consecutive episodes of size l, sharing the l -1 first events is called a block. Candidate episodes are identified by creating all possible combinations of two episodes of the same block [START_REF] Verkalo | Discovery of frequent episodes in event sequences[END_REF].

The collection of frequent episodes is a list featuring episodes that most frequently appear in the n sequences. The algorithm computes the collection of frequent episodes that make up the sequence s from the class of parallel episodes. The search is done by levels of length l on episodes α in the episode class following the sub-episode relationship. The search begins with the most general episode (one containing only one event). At each level, the algorithm calculates the collection of episodes and their frequencies.

The frequency of an episode is defined as the ration between the number of the window where the episode α appear entirely and the total window number. With a sequence of events s, a window of dimension win, the frequency of an episode α in s is:

f r(α, s, win) = |{w ∈ W (s, win)|α ∈ w}| |W (s, win)| (9)
where s is the number of windows is the α episode is fully included, W is the set of all windows w on s. For a defined threshold frequency thres_fr, α is frequent if f r(α, s, win) ≥ min_f r.

A rule for predicting frequent episodes that respect s, win, et min_f r is denoted F (s, win, min_f r). A prediction rule is a β ⇒ γ expression where β and γ are episodes such as β ≤ γ. The ratio f r(γ,s,win) f r(β,s,win) is called the confidence of the episode prediction rule. This confidence can be interpreted as a conditional probability that the entire γ is present in the window, with a given β.

Comparaison of LAE and the standard algorithm

In this section, we compare the algorithm LAE presented in section 3 to the standard algorithm of [START_REF] Verkalo | Discovery of frequent episodes in event sequences[END_REF] presented in section 5.1. The first step of the standard algorithm is to find all candidate episodes of the n sequences. Once all the episodes are found, they are stored and classed alphabetically in order to start creating all possible combinations of two episodes of the same block. Once the candidate episodes are generated, they are stored in the F (s, win, min_f r) episode collection. The prediction and confidence rules are then deduced from the previous result. The LAE algorithm proceeds quite differently: it begins by looking for events that are not in common in the n sequences and then it removes them from the sequences to keep only the common events. It determines the maximum occurrence number of each event, which allows us to easily find the maximum frequency of each pair of events on the n sequences.

The difference with the standard algorithm is that it finds the episodes in the form of a constraint graph describing all the possible episodes. The standard algorithm makes it possible to obtain an episode table directly. The algorithms have their differences, however: with the F (s, win, min_f r) event collection, it is possible to predict or anticipate the behavior of the sequence in its integrity. The algorithm thus provides a set of episodes as well as their prediction rules. On the other hand, the LAE algorithm provides important information thanks to the consecutive test. It makes it possible to express the fact that there can be no intermediate events between the two events of a pair on the n sequences of inputs. An intermediate event is an event of a different type than those defined in the pair. The comparison with the Mannila algorithm will be pursued when we'll have a larger set of data. Then, real tests and further analysis could be done.

Toward the Ph.D work

The thesis is the continuation of the presented work. The Ph.D is a collaboration with the ACTIA company which developed a software named ACTI-DIAG destined to mechanics. It allows to dialogue with the set of calculators (airbag, injectors, ...) of a vehicle and to pick up all the available information inside the car. These information can then be exploited in conjunction with information from the manufacturer or independent actors to finally allow the mechanic to locate the component to be replaced. In the field of automotive diagnostics, the number of identical vehicles circulating makes it possible to capitalize a great amount of knowledge. This capitalization then makes it possible to repair a vehicle more quickly by taking advantage of the experience gained on similar vehicles having already encountered the same problem. Currently, the information collected on the vehicle are:

• The type of the vehicle.

• The set of calculators et their actual software version.

• The set of default codes (also called Data Trouble Code) on calculators with the occurrence time. • The context (values of parameters judged relevant).

• etc... The objectives of the thesis is to analyze (semi-) automatically this set of data to reduce the updating time of the diagnosis and repair methods by the manufacturer, but also to reduce the time of the diagnosis and repair stages when a vehicle is in a workshop. This reduction can be addressed by:

• An update of existant knowledges (the list of possible causes for a given symptom, the probability of each fault occurrences about a vehicle, ...). • A ranking of knowledges proposed to the mechanic (most probable hypothesis related to a symptom, the "incident" files having most often completed, ...). • The creation of a new incident file about the most recurrent fails in the network. • Optimize existing diagnosis procedures to minimize their journey time by taking into account the time required to complete each of their tests and the probability of occurrence of each fault. • The learning of new diagnosis procedures by using collected data (actual DTCs, the context of the raise of DTCs, ...) in order to distinguish different possible faults for a or several given symptoms.

ACTIA is setting a network of 50 workshops for collecting data from different diagnosis with the symptoms and causes labels. This operation will allow to do real tests and comparisons with standard algorithms. The deployment of this network is planned for the end of November. Until that time, the study is performed without the label information. In first step, for a given ECU (Electronic Control Unit), the goal is to find co-occurrence of DTCs on vehicles. As DTC are discrete data (for example, P318D is the DTC for a communication error for the Nissan Leaf model), learning algorithms like clustering with non-ordered discrete data [START_REF] Watve | Clustering non-ordered discrete data[END_REF] and k-Nearest Neighbor Searching in non-ordered discrete data space [START_REF] Kolbe | On k-nearest neighbor searching in non-ordered discret data spaces[END_REF] are investigated.

Conclusion and perspectives

The work presented in this paper is part of a learning process for generic automotive diagnostics. The objective is to automatically generate a model of the diagnosis procedures used by mechanics to diagnose a given fault. This is performed from a set of action sequences carried out during several diagnosis sessions referring to the fault. Since the traces recorded by the diagnostic tool of our case study are timeless, we have directed ourselves to the development of our own learning algorithm: Learning Algorithm for Episodes. This algorithm allows to learn temporal patterns where time information is captured only through the notion of precedence. The tests carried out made it possible to validate the first version of the algorithm presented in the paper. This work has reinforced the interest of automatic learning to exploit the feedback from the diagnosis sessions. Several perspectives have been identified. First, it is necessary to test the algorithm with a larger volume of input data and to evaluate its efficiency. Also, the consecutive test need consolidation. From an application point of view, the idea is to use the results of LAE to improve the training tools of garage mechanics. As the complexity of vehicles increases, paper is not suffisant anymore. Our industrial partner (AC-TIA) develops a Serious Game to facilitate the learning of new methods of diagnosis and repair. A Serious Game is a software that combines a serious intention with playful springs. Episodes provided by LAE can be turned into training scenarios for the Serious Game for the diagnosis and repair tasks. The problem of knowledge capitalization addressed in this work is a real challenge for ACTIA company. The collaboration with ACTIA on this research area is going in the context of a Ph.D.

Figure 1 :

 1 Figure 1: Example of a serial episode and a parallel episode

2 )s 2

 22 1), (B, 2), (C, 3), (D, 4), (A, 5) > (=< (A, 1), (B, 2), (C, 3), (D, 4), (E, 5), (A, 6) > (3) s 3 =< (B, 1), (C, 2), (A, 3), (D, 4), (A, 5), (E, 6) > (4)

  1), (B, 2), (C, 3), (D, 4), (A, 5) > (6) s 2 =< (A, 1), (B, 2), (C, 3), (D, 4), (E, 5), (A, 6) > (7) s 3 =< (B, 1), (C, 2), (A, 3), (D, 4), (A, 5), (E, 6) > (8)

Figure 2 :

 2 Figure 2: Constraints graph of the 3 test sequences

Figure 3 :

 3 Figure 3: Overview of the context

Figure 4 :

 4 Figure 4: Extract of an Actia's sequence

Table 1 :

 1 List Com of common events in s 1 , s 2 et s 3Table 1 represents all the common event types located in S.Here, the common event types are A, B, C and D. Since event E is not present in all the sequences, it is not present in the list Com.

	Common events are searched in the
	input sequences then, sequences of entries are modified
	to keep only the common events. In this context, events
	are different screens used during the diagnosis session.
	Experienced mechanics receive training on how to diagnose
	a vehicle. They know how to use the diagnosis tool and
	which screens to navigate. The events that are not in all
	input sequences are considered as beginner mechanic's
	actions.
	A B C D

Table 2 :

 2 s 1 , s 2 et s 3 with only common events

	Table

Table 4 :

 4 Maximal frequencies in the table F req of each event types in s 1 , s 2 , s 3In the table 4, we can see the maximal frequency of each common event type. In this example, this information equals 2 for the event A and 1 for the others event types. To calculate the pairs available in each sequence, all combinaisons of F req[i] with F req[j] must be performed.

	AB AC AD BC BD CD
	2	2	2	1	1	1

Table 5 :

 5

Set of pairs available in sequences s 1 , s 2 , s 3

This work has been supported by the company ACTIA, one of the market leaders for diagnostic tools in the automobile domain.

Acknowledgements

We would like to thank our contacts in the company ACTIA for their valuable comments and recommendations.


Tom Obry

Louise Travé-Massuyès

Audine Subias

Computer-aided Diagnosis via Hierarchical Density Based Clustering

When applying non-supervised clustering, the concepts discovered by the clustering algorithm hardly match business concepts. Hierarchical clustering then proves to be a useful tool to exhibit sets of clusters according to a hierarchy. Data can be analyzed in layers and the user has a full spectrum of clusterings to which he can give meaning. This paper presents a new hierarchical densitybased algorithm that advantageously works from compacted data. The algorithm is applied to the monitoring of a process benchmark, illustrating its value in identifying different types of situations, from normal to highly critical.

Introduction

In data-based diagnosis applications, it is often the case that huge amounts of data are available but the data is not labelled with the corresponding operating mode, normal or faulty. Clustering algorithms, known as non-supervised classification methods, can then be used to form clusters that supposedly gather data corresponding to the same operating mode.

Clustering is a Machine Learning technique used to group data points according to some similarity criterion. Given a set of data points, a clustering algorithm is used to classify each data point into a specific group. Data points that are in the same group have similar features, while data points in different groups have highly dissimilar features. Among well-known clustering algorithms, we can mention K-Means [START_REF] Macqueen | Some methods for classification and analysis of multivariate observations[END_REF], PAM [START_REF] Kaufman | Clustering by means of medoids[END_REF], K-Modes [START_REF] Huang | Extensions to the k-means algorithm for clustering large data sets with categorical values[END_REF], DBSCAN [START_REF] Ester | A density-based algorithm for discovering clusters in large spatial databases with noise[END_REF].

Numerous validity indexes have been proposed to evaluate clusterings [START_REF] Halkidi | On clustering validation techniques[END_REF]. These are generally based on two fundamental concepts :

• compactness, the members of each cluster should be as close to each other as possible. A common measure of compactness is the variance, which should be minimized.

• separation, the clusters themselves should be widely spaced.

Nevertheless, one must admit that the concepts discovered by even the most scored clusterings hardly match business concepts [START_REF] Arbelaitz | An extensive comparative study of cluster validity indices[END_REF] [START_REF] Liu | Understanding of internal clustering validation measures[END_REF]. One of the reasons is that data bases are often incomplete in the sense that they do not include the data about all the influencial attributes. In particular, business concepts are highly sensitive to environmental parameters that fall outside the scope of the considered business domain and that are not recorded, for instance stock exchange. In addition, the clusters corresponding to business concepts may be quite "close" in the data space and the only way to capture them would be to guess the right number of clusters to initialize correctly the clustering algorithm. This is obviously quite hard. Hierarchical clustering then proves to be a useful tool because it exhibits sets of clusters according to a hierarchy and it modulates the number of clusters. Data can then be analyzed in layers, with a different number of clusters at each level, and the user has a full spectrum of clusterings to which he can give meaning. Hierarchical clustering identifies the clusters present in a dataset according to a hierarchy [START_REF] Stephen | Hierarchical clustering schemes[END_REF][9] [START_REF] Defays | An efficient algorithm for a complete link method[END_REF]. There are two strategies to form clusters, the agglomerative ("bottom up") strategy where each observation starts in its own cluster and pairs of clusters are merged as one moves up the hierarchy. The divise method ("top down") where all observations start in one cluster and splits are performed recursively as one moves down the hierarchy. The results of hierarchical clustering are usually presented in a dendrogram. A dendrogram is a tree diagram frequently used to illustrate the arrangement of the clusters. In order to decide which clusters should be combined or where a cluster should be split, a measure of dissimilarity between sets of observations is required. In most methods of hierarchical clustering, splits or merges of clusters are achieved by use of an appropriate metric like euclidean, manhattan or maximum distance.

Few algorithms propose a density-based hierarchical clustering approach like ↵-unchaining single linkage [START_REF] Martínez-Pérez | A density-sensitive hierarchical clustering method[END_REF] or HDBSCAN [START_REF] Campello | Density-based clustering based on hierarchical density estimates[END_REF]. In this paper, we present a new hierarchical clustering algorithm, named HDyclee, based on density that advantageously works from compacted data in the form hypercubes. This contribution is an extension of the clustering algorithm DyClee [START_REF] Barbosa | A novel algorithm for dynamic clustering: Properties and performance[END_REF], [START_REF] Barbosa | A data-based approach for dynamic classification of functional scenarios oriented to industrial process plants[END_REF], [START_REF] Barbosa | A data-based dynamic classification technique: A two-stage density approach[END_REF]. The purpose of this work is to generate a flat partition of clusters with a hypercube's density level higher or equal to a threshold and to be able to visualize all existant clusters in the dataset with a dendogram by varying the density of the hypercubes present in a group. The value of the algorithm in a diagnosis context is illustrated with the monitoring of a Continuous Stirred Tank Heater benchmark, for which it allows the user to identify different types of situations, from normal to highly critical. This paper is organized as follows. In section 2 the Dy-Clee algorithm is presented. In section 3 the concepts and principles underlying Dyclee, like the definition of micro clusters µC, dynamic clusters and the KD-Tree structure, are explained. In the section 4, the hierarchical clustering based-density algorithm is presented. Tests and results are detailed in section 5. The conclusion and perspective for future work end this paper in section 6. 12 Dyclee: a dynamic clustering algorithm DyClee is a dynamic clustering algorithm which is able to deal with large amounts of data arriving at fast rates by adopting a two stages strategy similar to [START_REF] Charu C Aggarwal | A framework for clustering evolving data streams[END_REF], [START_REF] Kranen | The ClusTree: indexing micro-clusters for any stream mining[END_REF], [START_REF] Cao | Density-based clustering over an evolving data stream with noise[END_REF]. The first stage is a fast scale distance-based algorithm that collects, pre-processes and compresses data samples to form so-called micro-clusters (µ-clusters). It operates at the rate of the data stream and creates µ-clusters putting together data samples that are close, in the sense of a given distance, to each other. µ-clusters are stored in the form of summarized representations including statistical and temporal information. The second stage is a slower scale density-based algorithm that groups the µ-clusters into actual clusters that can be interpreted semantically as classes. It takes place once each t slow seconds and analyses the distribution of µ-clusters. The density of a µ-cluster is considered as low, medium or high and is used to create the final clusters by a density based approach, i.e. dense µ-clusters that are close enough (connected) are said to belong to the same cluster. Similarly to [START_REF] Chen | Density-based clustering for real-time stream data[END_REF], a cluster is defined as the group of connected µ-clusters where every inside µ-cluster presents high density and every outside µ-cluster exhibits either medium or low density. The above dense µ-cluster structure allows the algorithm to create clusters of non convex shapes even in high dimensional spaces and it has proved outliers rejection capabilities in evolving environments [START_REF] Cao | Density-based clustering over an evolving data stream with noise[END_REF]. In addition, µclusters of similar densities can form clusters of any shape and any size. In DyClee, both stages work on-line, but operate at different time scales. This multi-density feature allows the detection of novelty behavior in its early stages when only a few objects giving evidence of this evolution are present. Figure 1 gives the global description of DyClee.

Main principles of DyClee

All the principles explained in this section are from the core algorithm.

Notion of micro-clusters µC

Considering a d-dimensional object X = [x 1 , ..., x d ] marked with a timestamp t X and qualified by d features, a µ-cluster gathers a group of data samples close in all dimensions and whose information is summarized in a characteristic feature vector (CF). For a µ-cluster µC k , CF has the following form:

CF k = (n k , LS k , SS k , t lk , t sk , D k , Class k ) . (1)
where n k 2 < is the number of objects in the µ-cluster k, LS k 2 < d is the vector containing the linear sum of each feature over the n k objects, SS k 2 < d is the square sum of features over the n k objects, t lk 2 < is the time when the last object was assigned to that µ-cluster, t sk 2 < is the time when the µ-cluster was created, D k is the µ-cluster density and Class k is the µ-cluster label if known. Using LS k , SS k and n k the variance of the group of objects assigned to the µ-cluster k can be calculated.

The µ-cluster is shaped as a d-dimensional box since the L1-norm is used as distance measure. The distance between an object X = [x 1 , ..., x d ] T and a µ-cluster µC k , named as dis(X, µC k ), is calculated as the sum of the distances between the µC

k vector center c k = [c 1 k , . . . , c d k ]
T and the object value for each feature as shown in equation ( 2):

dis(X, µC k ) = L 1 (X, c k ) = d X i=1 x i c i k . ( 2
)
The data is normalized according to the data context, i.e. the feature range [min i , max i ] of each feature i, i = 1, ..., d. If no context is available in advance, it may be established online. The size of the hyperboxes S i along each dimension i is set as a fraction of the corresponding feature range. The hyperbox size per feature is hence found according to [START_REF] Huang | Extensions to the k-means algorithm for clustering large data sets with categorical values[END_REF], where i is a user constant parameter in the interval (0, 1), establishing the fraction:

S i = i |max i min i |, 8i = 1, . . . , d. (3)
Whenever an object X arrives, the algorithm searches for the closest µ-cluster. Once found, a maximal distance criterion is evaluated to decide whether or not the object fits inside the µ-cluster hyper-box. If the fitting is sufficient the µ-cluster feature vector is updated using the object information; if not, a new µ-cluster is created with the object information using its time-stamp as cluster time of creation.

The density of a µ-cluster µC k is calculated using the current number of objects n k and the current hyper-volume of the bounding box

V k = Q d i=1 S i , as shown in (4): D k = n k V k . ( 4
)
Let µC k↵ and µC k be two µ-clusters, then µC k↵ and µC k are said to be directly connected if their hyper-boxes overlap in all but ' dimensions, where ' is an integer. The parameter ', fixed by the user, establishes the feature selectivity. A µ-cluster µC k1 is said to be connected to µC kn if there exists a chain of µ-clusters {µC k1 , µC

k2 , • • • , µC kn } such that µC ki is directly connected to µC ki+1 for i = 1, 2, • • • , n 1.
A set of connected µ-clusters is said to be a group.

Dynamic clusters

Dyclee is a dynamic clustering algorithm, which means that not only the parameters but the classifier structure changes according to input data in an automatic way. It achieves several cluster operations like creation, elimination, drift, merge, and split. For instance, a cluster is splitted into two or more clusters if, with the arrival of new data, high density regions can be distinguished inside the cluster. In that scenario, dense regions are separated by low density regions, making the cluster no longer homogeneous. Even more, the cluster center could be situated in a low density region, loosing its interpretability as prototype of the elements in the cluster. Splitting the cluster creates smaller homogeneous clusters, completely representative of the belonging samples. An illustrative example of this phenomenon is shown in Figure 2.

Finding groups of µ-clusters: the KD-Tree structure

To find groups of connected µ-clusters, a KD-Tree [START_REF] Maneewongvatana | It's okay to be skinny, if your friends are fat[END_REF] is used. A KD-Tree is a binary tree, each of whose nodes represents an axis-aligned hyperrectangle. Each node specifies an axis and splits the set of points based on whether their coordinate along that axis is greater than or less than a particular value. The tree is queried to return only neighbors who are at a maximum distance from a point. A µC j is the neighbor of the µC k if the condition in equation ( 5) is respected :

L 1 = d max i=1 |x i k c i j | r . ( 5
)
where d is the number of dimension, c i j the center of the µC j at the dimension i and r the maximal distance from the µC k . In this context, r is set to i .

A new hierarchical clustering density-based algorithm

In this section, a new hierarchical clustering density-based algorithm is presented. Inputs are the connections between µ-clusters from the KD-Tree and the output is a flat partition of clusters where all µ-clusters that are in clusters have a minimum density level guaranteed.

Representation of the µ-clusters connections

The representation of the connections between all µ-clusters can be visualized by a weighted Graph. A weighted Graph G = (N , E, W) is a triplet where N is a set of nodes. A node n i corresponds to the µ-cluster µC i . E is the set of edges with e ij the edge between the node n i and n j , which are unordered pairs of elements of G. Finally, W is a set of weights on E with w ij defined in the equation 6.

w i,j = min(D i , D j ). ( 6
)
where the density of a µ-cluster D i is defined in the equation 4.

An edge e ij between µ-clusters µC i and µC j means those are directly connected in the sense defined in the section 3.1. If two µ-clusters are not directly connected, there is no edge between them which leads to a Graph that is not full. The Graph of µ-cluster's connections is built according to the algorithm 1. Neighbors are searched for each µ-cluster with respect to the equation 5 (lines 3 to 6). The function Search_neighbors() is detailed in the algorithm 2. For each µC k , the distance L 1 defined in equation 5 is applied where r = i , x i k the value of µC k at the i th dimension and c i j the center of the µC j at the dimension i. The variable Neigh-bors_of_k contains all the neighbors of µC

k . An edge e kj is added to the Graph G for each neighbor µC j of the µcluster studied µC k and the weight w kj is calculated with the equation defined above (lines 7 to 11). The weights are calculated for every edge following the equation 6 : (7)

The Figure 5 shows the weighted Graph which represents the connections between all µ-clusters in the dataset.

Representation of the hierarchy of clusters

The objective of the algorithm proposed in this paper is to represent and visualize all the possible clusters at different levels of density in a dataset. In contrast to most known hierarchical clustering algorithm, links that relies level l n and l n + 1 in our new algorithm's dendogram are not based on the distance between objects but on their densities. Furthermore, our proposal is not based on the objects but on µ-clusters that contains the objects to decrease the complexity of calculation. At the root of the tree, there is one cluster composed by all µ-clusters. Each cut in the tree corresponds to a density threshold, i.e each cluster formed below this cut level is composed by µ-clusters that have a density higher to the density threshold. At the bottom of the tree, there is one cluster for each µ-cluster. So root's density level is 0 and the last cut on the tree is max(w 2 W). Like [START_REF] Campello | Density-based clustering based on hierarchical density estimates[END_REF], a variable " is user parameter evolving in the interval " 2 [0, max(w 2 W)]. For each iteration of ", all the weights w 2 W are checked. if w ij "

, the edge e ij is removed. The vertical axis of the dendogram is 1/D to keep an ascending direction. Figure 6 represents the case when " = 3. We can observe two clusters composed by µ-clusters that have their densities strictly higher to three and one µ-cluster alone. This method allows to detect a split of clusters case and to isolate the least denses µ-clusters as early as possible. Once all values in the interval of " studied and all edges removed, the dendogram can be generated (see Figure 7).

Extracting a flat partition of clusters

The dendogram generated, a partition of clusters corresponding to a specific density level can be extracted from the Graph of µ-cluster's connections with a density threshold according to the algorithm 2. To find clusters, the edges that have a weight strictly lower than the density threshold are removed (lines 1 to 7). The remaining edges in the Graph have their weights higher or equal to the density threshold hence µ-clusters forming the clusters are guaranteed to have their densities higher or equal to the threshold. Remaining groups are searched in the Graph. If the size of a group is equal to 1, i.e the group have not connection with other µ-clusters, it is considered as noise. Else, the group is recognized as a cluster (lines 10 to 16). In this example, the density threshold is fixed to " = 8. Every group of µ-clusters that is below the density threshold is considered as a cluster. If a µ-cluster does not have connection to the other µ-clusters, then it is considered as noise. The Figure 10 illustrates how to visualize the clusters that have a density strictly above the threshold. Once the density threshold known, all the edges that have their weights strictly lower are removed to left only groups of µ-cluster to form final clusters. Figure 11 shows the final clusters. The first includes µC 1 and µC 2 , the second includes µC 4 , µC 5 and the last contains µC 6 , µC 7 . µC 2 and µC 8 are considered as noise because they do not have any connections with the other µ-clusters. for all w(i, j) in G do if w(i, j) " then

G.removed(e(i,j))

5:
end if

6:
end for 7: end for 8: Groups = Find_Groups() 9: k = 0 10: for all g in Groups do HDyClee is tested on a benchmark similar to the well known Continuous Stirred Tank Heater (CSTH) of [START_REF] Thornhill | A continuous stirred tank heater simulation model with applications[END_REF]. The CSTH is a stirred tank in which hot and cold water are mixed and further heated using steam. The final mix is drained using a long pipe [START_REF] Barbosa | A data-based approach for dynamic classification of functional scenarios oriented to industrial process plants[END_REF]. Figure 12 For this experiment, the radius of research for a µcluster's neighbors r is set to 0.06. The parameter ' that defines the number of dimensions that must overlap so that two µ-clusters are considered directly connected is set to 0. That means two µ-clusters µC k↵ and µC k are directly connected if their hyperboxes overlap in all dimensions. The parameter " is set to 0 for results shown in Figure 16 that correspond to the root of the dendogram. For the following graphs, the x-axis is the flow CW flow normalized and the y-axis is the tank temperature T ank temperature normalized. The tank level T ank level is not plotted. Figure 14 shows the graph of connections between µ-clusters when " = 0. Each red square represents a µ-cluster. Micro-clusters that are not connected to the others are considered as noise. The µcluster µC 293 is connected to µC 222 and µC 232 , meaning there are directly connected. µC 222 is connected to µC 34 and µC 232 is connected with µC 41 and so on. This chain of µ-clusters forms a cluster. Figure 15 shows a generalized dendogram representing the hierarchy of behaviors found in the dataset. It is possible to visualize the most frequent behaviors of the system, in our case the normal behavior and fault l 3 . This is reproted in Figure 17. For this purpose, the density threshold is set to " = 6 by using the dendogram. At this density, the clusters corresponding to other behaviors are considered as noise because their maximal densities are less than 6. To visualize only the nominal behavior, the dendogram must be cut at the density threshold " = 7 because above this value, the other clusters have no µ-clusters that are connected to each other. This is shown on Figure 19. Figure 18 illustrates the graph of µ-cluster connections after deletion of the edges with a weight w ij 7

. Some µ-clusters which were part of the biggest cluster are now considered as noise. Indeed, edges that connected them to other µ-clusters were less than the density threshold.

Conclusion and perspectives

The work presented in this paper proposes a new hierarchical density-based algorithm named HDyclee. The purpose of this algorithm is to extract a hierarchy of clusters that are guaranteed to have a level of density at each layer. Branches in the dendogram do not represent distance between objects but minimum density difference. This approach allows one to identify clusters with the poorest densities and then walk up the hierarchy for higher densities. The algorithm is detailed and tested on a well known monitoring benchmark. HDyClee is able to detect all the behaviors of the process and the user can explore more or less frequent behaviors by cutting the dendogram at different densities.

Next step is to develop experimentations in order to compare this new algorithm with other density-based algorithms. Then the comparative study will include hierarchical and distance-based clustering methods [START_REF] Rodríguez Ramos | An approach to fault diagnosis using fuzzy clustering techniques[END_REF], [START_REF] Lei | New clustering algorithm-based fault diagnosis using compensation distance evaluation technique[END_REF]. Several perspectives have been identified for HDyClee, which follow from DyClee properties. In particular, DyClee has a forgetting function that allows to forget µ-clusters which do not receive any data or are not significant (not denses µclusters). This function will be included in HDyClee, which will allow us to produce a dynamic dendogram and then to visualize the evolution of the different behaviors of a system.

Figure 1 :

 1 Figure 1: Global description DyClee.

Figure 2 :

 2 Figure 2: Splitting of cluster A.

Algorithm 1 2 1 :

 121 Build the Graph of µ-cluster's connections Require: KD-Tree 1: G = Graph() 2: Connection = () 3: for k = 1 to Number of µ-clusters do Research of a µ-cluster's neighbors Require: KD-Tree, µ-cluster µC k for j = 1 to Number of µ-clusters do end for 7: return Neighbors_of_k Let us consider five µ-clusters µC 1 , µC 2 , µC 3 , µC 4 and µC 5 with their densities D 1 = 14, D 2 = 12, D 3 = 3, D 4 = 9, D 5 = 13. Figure 3 shows those five µ-clusters.

Figure 3 :

 3 Figure 3: Example of five µ-clusters.

Figure 4 :

 4 Figure 4: Search of neighbors for the µC 1 .Once neighbors of µC 1 are found, the neighbors of the other µ-clusters are searched. The result is shown in table1:

w

  12 = min(D 1 , D 2 ) = min(14, 12) = 12, w 23 = min(D 2 , D 3 ) = min(12, 3) = 3, w 34 = min(D 3 , D 4 ) = min(3, 9) = 3 w 45 = min(D 4 , D 5 ) = min(9, 13) = 9.

Figure 5 :

 5 Figure 5: Example 1: Graph of connections between µclusters.

Figure 6 :

 6 Figure 6: Example 1: Graph of connections between µclusters after removed all weight's edges w ij < 3.

Figure 7 :

 7 Figure 7: Dendogram of the example 1.

Figure 8 :

 8 Figure 8: Example 2 : Graph of connections between eight µ-clusters.

Figure 9 :

 9 Figure 9: Dendogram of the figure 8.

Figure 10 :

 10 Figure 10: Visualization of clusters in the dendogram that satisfies the " = 8.

Figure 11 :

 11 Figure 11: Clusters formed with a density guaranteed strictly above 8.

  shows the structure of the CSTH. Process inputs are set-points for the cold water, hot water and steam valves. Process outputs are hot and cold water flow, tank level and temperature. Process inputs and outputs represent electronic signals in the range 4-20 mA. The test is done using three output variables: cold water flow CW flow , tank level T ank level , and temperature of the water in the tank T ank temperature , in the operation mode OP 1. In this mode, these variables are regulated at the values provided in

1 1 fixed l 2 Al 1 l 2 3 :

 11223 , and T ank temperature are shown in Figure 13 and their recorded values across time constitute the data set for our hierarchical clustering experiment. Sudden changes in the value of regulated variables are indicative of the occurence of some fault or of some fault being fixed. The dataset was generated by simulation. Event Description l Evolving leak starts. Hole diameter goes from 1 to 3, 5mm in 1500 seconds l Leak second evolving leak starts. The second hole goes from 0 to 1mm in 1500 seconds Leaks Description of faults on the CSTH system for the operation mode OP 1.

Figure 13 :

 13 Figure 13: Process measurements for multiple fault scenario.

Figure 14 :

 14 Figure 14: Graph of µ-clusters connections for the study case when " = 0.

Figure 15 :

 15 Figure 15: The dendogram with the nominal behavior and the occurrence of faults.

Figure 16 :

 16 Figure 16: Clusters found by HDyClee algorithm: the nominal behavior (green cluster) and abnormal situations.

Figure 17 :

 17 Figure 17: The nominal behavior (blue) and the most frequent fault (green).

Figure 18 :

 18 Figure 18: Graph of µ-clusters connections for the study case after removing the weight's edges w ij < 7.

Figure 19 :

 19 Figure 19: Visualization of the nominal behavior after having cut the dendogram at level " =7.

table 1 :

 1

	µ-clusters Neighbors
	µC 1	µC 2
	µC 2	µC 1 , µC 3
	µC 3	µC 2 , µC 4
	µC 4	µC 3 , µC 5
	µC 5	µC 4

Table 1 :

 1 Neighbors for each µ-clusters

table 2 .

 2 The process undergoes several Figure 12: The continuous stirred tank heater.faults and several repairs that are reported in table3.

	Variable	OP 1
	CW flow (mA)	11.89
	T ank	level (mA)	12.00
	T ank

temperature (mA) 10.50

Table 2 :

 2 Nominal values for the test on CSTH.

	The measurements of CW

flow , T ank level

This work is performed in the framework of a CIFRE Project supported by ACTIA.


J Margueron

N Van Giai

EQUATION OF STATE IN THE INNER CRUST OF NEUTRON STARS: DISCUSSION OF THE UNBOUND NEUTRONS STATES

Keywords: Equation of state, self-consistent mean-field model, pairing correlations

In this paper, we calculate the stable Wigner-Seitz (W-S) cells in the inner crust of neutron stars and we discuss the nuclear shell effects. A distinction is done between the shell effects due to the bound states and those induced by the unbound states, which are shown to be spurious. We then estimate the effects of the spurious shells on the total energy and decompose it into a smooth and a residual part. We propose a correction to the Hartree-Fock binding energy in Wigner-Seitz cell (HF-WS).

Introduction

From the very first nuclear models based on extrapolations of a mass formula, 1 several models have been developed to provide the equation of state in the inner crust, like the BBP 2 and BPS 3 semi-classical models or fully quantal self-consistent mean-field models. [4][5][6][7]9,10 The latter models are based on the Wigner-Seitz (W-S) approximation and the continuum states are discretized in a box (spherical or not). In semi-classical models, a clear limitation is coming from the absence of nuclear shell effects which is known to play an important role in determining the composition of the ground state prior to the neutron drip, 4 in the inner crust 6 as well as in the transition zone between the crust and the core. 11 Since the seminal work of Negele & Vautherin, 4 many efforts have been invested in the developpement of self-consistent mean-field modelization of the inner crust, by introducing deformation of the nuclear cluster in the high density part of the crust, 5,10 or pairing correlations, [7][8][9] or more recently by improving of the lattice description within the band theory. 12 Band theory takes into account the proper symmetries of the system, but the equations are numerically very complicated to solve and have not yet been solved in a self-consistent framework. Instead, it has been coupled to a selfconsistent mean-field model at the W-S approximation. 13 It has then been shown that the W-S approximation is justified if the temperature is larger than about Table 1. Comparison of the proton number Z and the radius of the W-S cell R W S obtained in different calculations. See the text for more explanations. 100 keV, or if the quantity of interest averages the density of states on a typical scale of about 100 keV around the Fermi surface. The typical scale which has been found, 100 keV, is in fact related to the average inter-distance energy between the unbound states which are discretized in the W-S box. In this paper, we clarify the effect of the discretization of the unbound states on the ground state energy as well as the role of pairing correlations. Consequences for the equation of state are presented.

Nzone

k F Z R W S ( in fm) fm -1 N-V 4 HF-BCS a HF HFB N-V

Equation of state in the inner crust

The inner crust matter is divided in 11 zones shown in Tab. 1, as in Ref. 4 (notice that the denser zone has been removed since it is in the deformed pasta region). Each W-S cell is supposed to contain in its center a spherical neutron-rich nucleus surrounded by unbound neutrons and immersed in a relativistic electron gas uniformly distributed inside the cell. For a given baryonic density ρ B = A/V , neutron number N and proton number Z, the total energy of a W-S cell, E tot has contributions coming from the rest mass of the particles, E m (N, Z) = Z(m p + m e ) + N m n , from the nuclear components (including Coulomb interaction between protons) E n (ρ B , N, Z), from the lattice energy,

E l (ρ B , N, Z, {ρ p }) = 1 2 d 3 r 1 d 3 r 2 ρ e (r 1 ) e 2 |r 12 | (ρ e (r 2 ) -ρ p (r 2 )) , (1)
which is induced by the difference between the uniform electron density and proton density, and finaly from the kinetic energy T e of the relativistic electron gas. 4 The Coulomb exchange energy as well as the screening correction for the electrons gas has been neglected. The self-consistent Hartree-Fock-Bogoliubov (HFB) approach has been presented in detail in Ref. 14 It has been extended to describe the dense part of the inner crust at finite temperature in Ref. 15 successfully applied to calculate the specific heat of the whole part of the inner crust. 16 The effective N-N interaction should provide a reasonable description of both the nuclear clusters and the neutron gas. We use the Skyrme force SLy4 17 which was adjusted to describe properly the mean field properties of neutron-rich nuclei and infinite neutron matter. The choice of the pairing force is more problematic since at present it is not yet clear what is the strength of pairing correlations in low density neutron matter. However, in the present paper we have chosen to adjust the effective pairing interaction on the results of the pairing gap obtained with the bare interaction. We then have a maximum gap in neutron matter of about 3 MeV located at k F = 0.85 fm -1 . The minimization is performed in two steps: first, at a given density ρ B and for a given mass number A, we look for the cell which satisfies the beta-stability criterion, µ p + µ e = µ n . The neutron and proton chemical potentials are provided by the HFB calculation while the electron chemical potential is deduced from µ e = dE tot /dN e . Note that, due to shell effects we may find several cells which satisfy this criterion. We then choose the one which has a minimal energy. In the second step, we minimize the total energy B tot with respect to the mass number A. It should be noted that, in the second step all the cells must be calculated at the same density while the minimum energy is searched over the variable A. As the volume is modified by unit steps (the step of the mesh is 0.2 fm), the mass number A can take only discrete values so that A = ρ B V is always satisfied. The step between two values of A is not constant and vary as the square of the WS radius R W S . We have approximatively dA ∼ ρ B 4πR 2 W S dR W S . For N zone =4, we show in Fig. 1 a typical example of step 2 in the search for the W-S cell which minimizes the energy. On the l.h.s. is represented the total binding energy E tot /A of the W-S cells obtained in Hartree-Fock (HF) mean-field model (without pairing). The position of the lowest minimum is indicated as well as two other local minima. In this case, there is a difference of 5 keV per nucleon between the two first minimum. On the r.h.s., we plot the results of HFB calculations which show the effects of the pairing correlations: the first minimum is deeper (about 34 keV per nucleons). On the lower part of the graph, we show the partial contributions of the rest mass energy E m /A -m n , the nuclear binding energy E n /A, the electron kinetic energy T e /A and lattice energie E l /A. One sees that the inner crust is a typical example of frustrated system: when the nuclear energy is minimum, the electron energy is maximum and vice-versa. The minimum in energy is a compromise between these two opposite energies. The consequence is that it gives rise to numerous local minima which are not very far in energy from the absolute minimum.

The equation of state obtained at the HF and HFB approximations is given in Table 1. In total, 19714 cells configurations has been calculated. The pairing correlations do not strongly modify the HF results, it may slightly decrease the proton number and the WS radius. We also show previous calculations of equations of state. We have obtained proton numbers Z different from those of Refs. 4,9 The results of Ref. 4 have been obtained in a density matrix expansion framework interpolating microscopic calculations for homogeneous matter. The difference may be due to the fact that the Skyrme interaction SLy4 is based on more recent microscopic parametrizations. It is difficult to compare with the results of Ref. 9 since those calculations have not been done at the same densities and they have essentially been carried out in the dense part of the crust. We then show an interpolation of their results. The comparison shows a similar trend in the dense cells.

Discussion of the shell effects

In ordinary nuclei, the bound states play an essential role in determining the ground state properties and the role of the continuum states is usually negligeable. The coupling between bound states and continuum states can be more important for nuclei near the drip lines, as the Fermi energy is increasing. However, even in these most extreme situations continuum states are only virtually populated and the occupation numbers associated to continuum states are small compared to those of the bound states. In the inner crust of neutron stars, the situation is completely different: the Fermi energy is positive and most of the mass is due to the unbound neutrons while bound states account only for a small fraction. Thus, the continuum states contribute in an essential way to physical properties of the crust like its density, total energy or specific heat. The presence of bound states, or nuclear clusters, is not negligible but it induces only corrections to the ground-state properties. 16 It is then very important to have an accurate model of the unbound states.

It should also be noted that, as the continuum states are embedded in a Coulomb lattice, their asymptotic behavior is different from that occurring in nuclei. The continuum states are indeed periodically distorted within a distance equal to the lattice step (20-100 fm). The proper theory for the description of the continuum states is the band theory which requires to solve the Schrödinger equation, 12,13 (h 0 + h k )u α,k (r) = ǫ α,k u α,k (r), where h 0 is the usual single-particle Hamiltonian of mean-field models. The second term h k is induced by the Floquet-Bloch boundary conditions for the single-particle wave function ϕ α,k : ϕ α,k (r) = u α,k (r)e ik•r . The spherical W-S approximation is obtained by setting k = 0 and replacing the W-S elementary polyhedron by a sphere. Then, the remaining index α is discrete, like in usual mean-field models, and runs over the eigenstates of the Hamiltonian h 0 . The main difference between mean-field models and band theory is indeed due to the term h k . In band theory, this term introduces a new index k which is a continuous variable ranging between 0 and π/a fm -1 where a is the lattice step. As a consequence, the density of unbound states is continuous (notice that it may also have structures and gaps). It is then clear that the discrete distribution of unbound states present in usual mean-field models is a consequence of the dropping of h k in the W-S approximation.

We can characterize the distribution of unbound states by the average distance between the unbound states energies, which can be related to the radius of the W-S cell as ∆ǫ ∼ /2mR 2 W S : the smaller the W-S radius, the larger the average distance ∆ǫ. One could notice that the W-S radius is a function of the density and it is fixed by the equation of state, as shown in Table 1. At the W-S approximation, the distribution of unbound states is then also fixed by the density. The non-continuous distribution of unbound states could be interpreted as a spurious effect, induced by the W-S approximation, since in band theory calculations, the density of unbound states varies in a continuous way. 13 One can then distinguish between the physical shell effects due to the bound states from the spurious shell effects due to the unbound states and the W-S approximation.

Spurious shell effects

In order to estimate the spurious shell effects due to the unbound neutrons, we simulate an homogeneous gas of neutron matter in W-S cell by removing all the proton states. In homogeneous matter, the energy density should be independent of the volume, then of the W-S radius, and it should depend only on the density. In our simulation of homogeneous matter, we then vary the W-S radius from 10 to 50 fm, at constant neutron density. We call B W S-hom. (ρ unb. , R W S ) the binding energy per particle obtained for a given W-S radius R W S and for a given density of unbound neutrons ρ unb. . We show the binding energy versus the W-S radius in Fig. 2(a) for several densities corresponding to N zone =6, 7 and 8. It is clearly seen that the binding energy is not independent of the W-S radius, except for large W-S radii where the binding energy converge to the energy of homogeneous matter, B hom. (ρ unb. ). The binding energy of homogeneous matter B hom. (ρ unb. ) is calculated

10 20 30 40 50 R WS [fm] 0,5 1 1,5 2 2,5 B
WS-hom. (ρ unb. ,R WS ) [MeV] 10 20 30 40 50 R WS [fm] -0,2 0 0,2 0,4 0,6 0,8 B hom. (ρ unb.
) -B WS-hom. (ρ unb.

,R WS )

10 20 30 40 50 R WS [fm] -0,2 0 0,2 0,4 0,6 0,8 B hom. (ρ unb.
) -B

WS-hom. (ρ unb. ,R WS ) -f(ρ unb. ,R WS ) (a) (b) (c)
N zone = 6

N zone = 7

N zone = 8 As expected, this quantity is converging to zero for large W-S radii, but for small radii, around 20 fm for instance, the difference can be as large as 300 keV per nucleons.

From the dots represented on Fig. 2(b), we obtain an universal function which fits the smooth radial dependence of the difference B hom. (ρ unb. ) -B W S-hom. (ρ unb. , R W S ). As this effect is due to the discretization of the unbound states, it is proportional to the average distance ∆ǫ between the box states. This difference should then decrease with the W-S radius as R -2 W S . It should also be proportional to the number of unbound neutrons, while the density dependence is not known. We then represent it by a power law ρ α unb. where the power α is adjusted to reproduce the dots presented in Fig. 2(b). In this figure we show the fitted function,

f (ρ unb. , R W S ) = 89.05(ρ unb. /ρ 0 ) 0.1425 R -2 W S . (2)
It interpolates the smooth part of the shell effects on the binding energy. There is, however, a residual difference which cannot be fitted. Finally, to estimate the residual difference between B hom. (ρ unb. ) -B W S-hom. (ρ unb. , R W S ) and the adjusted function f (ρ unb. , R W S ), shown on the r.h.s of Fig. 2(c). The average residual fluctuations are now of the order of 50 keV.

Modified HF binding energy in W-S cells (HF-WS)

In Sec. 4, it has been shown that the fluctuations in the binding energy induced by the spurious shell effect can be decomposed into a smooth and a residual term and we have obtained a universal function for the smooth term. We then propose a systematic correction to the HF binding energy in the W-S cell:

B HF -W S (ρ B , ρ unb. , R W S ) ≡ B W S-inhom. (ρ B , R W S ) + B hom. (ρ unb. ) -B W S-hom. (ρ unb. , R W S ) smooth = B W S-inhom. (ρ B , R W S ) + f (ρ unb. , R W S ) . (3)
In the calculation of the total energy, the nuclear binding energy E n /A should then be replaced by B HF -W S (ρ B , ρ unb. , R W S ). The minimum energy can be unambiguously identified if the difference in energy between the first and the second minimum is larger that the residual fluctuation, about 50 keV. In the HF calculation, the difference in energy between different cell configurations are of the order of few keV to few tens of keV. The residual fluctuation of the proposed correction is then too large. However, as seen above, the pairing correlations help in producing a deeper first minimum. It may also help in reducing the size of the residual fluctuations. In a future investigation, we plan to obtain an improved binding energy for W-S cells including pairing correlations. It must be remarked that the effects of the coupling between the protons and the unbound neutrons have been neglected in B HF -W S . This is a good approximation since the protons are deeply bound.

Conclusions

In the search for the W-S cells which minimize the energy, it is important to look not only at the lowest minimum energy, but also at the other minima. The differences between the first and other minima show how stable is the lowest minimum with respect, for instance, to thermal fluctuations or spurious shell effects. In HF calculations, the differences in total energy between the first minimum and the others are about few keV per nucleon (few MeV in total energy). It is shown that the spurious shell effects, induced by the W-S approximation, introduce a fluctuation in the calculation of the total energy up to about 300 keV per nucleon. Interpolating the smooth part of the spurious shell effects, we have then proposed a simple method to reduce the fluctuations down to about 50 keV per nucleon (HF-WS). The residual fluctuations are still too large, and the removal of the smooth part is not enough to improve the results obtained within the W-S approximation. We have also shown that the pairing correlations help in producing a deeper first minimum. We have then compared the equation of state obtained from HF and HFB binding energies. Several issues should still be addressed in future studies: spurious shell correction with pairing correlations as well as effects of the temperature. These two effects should give lower residual fluctuations than the one we obtained at the T = 0 HF approximation.

Fig. 1 .

 1 Fig.1. The binding energy ǫ/ρ B of the W-S cells in zone Nzone=4 as well as the partial contributions of the rest mass energy Em/Amn, the nucleons, the electrons and the lattice. HF and HFB results are shown on the left and right panels, respectively. See the text for more details.

Fig. 2 .

 2 Fig. 2. We represent several quantities versus the W-S radius: (a) the binding energy B W S-hom. (ρ unb. , R W S ) corresponding to Nzone=6, 7 and 8 (b) the difference B hom. (ρ unb. ) -B W S-hom. (ρ unb. , R W S ) (dots) and the function f (ρ unb. , R W S ) which fit the smooth component (solid lines), and (c) the difference between the dots presented in (b) and the function f (ρ unb. , R W S ).