Keywords: System Engineering, Dependability, Redundancy Policy, Phased Mission System, UML Class Diagram, Markov Chain

Systems engineering (SE) is a very promising approach to facilitate the development of complex systems. This explains why several SE processes have been already proposed. However, these proposals focus mainly on systems with faultless components. Integration of dependability concerns into SE processes must be supported by a suitable organization of the data which are dealt with during the system life-cycle. A metamodel which defines the concepts used during this cycle as well as the relations between these concepts is a way to rigorously describe this organization.

This article proposes such a meta-model developed for power production systems. These systems are phased mission systems composed of repairable and multi-state components; moreover, several redundancy policies shall be defined for each phase. This proposal is illustrated on a small example from a power plant. Last, the merit of this contribution to support the integration of dependability concerns is shown by proposing a method to build systematically, from the instance diagrams derived from the proposed meta-model, the Markov Chains which represent the dysfunctional dynamic behavior of a system.

I. INTRODUCTION

N the current socio-economical context, where costs, delays and dependability are crucial concerns while promising technological solutions frequently appear, engineering of critical systems is a complex issue that must be thought in the framework of systems engineering [START_REF] Leveson | Engineering a Safer World: Systems Thinking Applied to Safety[END_REF] and [START_REF] Haimes | Risk Modeling, Assessment, and Management[END_REF]. This approach permits in particular to avoid that incorrect industrial practices, which focus on counting or managing failures instead of preventing them, as detailed in [START_REF] Barnard | What is wrong with Reliability Engineering[END_REF], are introduced during the development of the system. The numerous activities of the systems engineering process (requirements analysis, functional analysis, risk management, dependability analysis, verification and validation, etc.) must be supported nevertheless by a suitable organization of the data they produce or consume. An efficient solution to describe this organization is to develop a meta-model, in the form of a P.-Y. Piriou and J.-M. Faure are with Ecole Normale Superieure of Cachan, France (e-mail: {pierre-yves.piriou; jean-marc.faure}@lurpa.enscachan.fr) G. Deleuze is with Electricite de France R&D, Clamart, France (e-mail: gilles.deleuze@edf.fr) UML/SysML class diagram for instance. A meta-model defines the concepts which are used during the engineering process as well as the relations between these concepts, then ensures data consistency, and facilitates the automatic construction of some models from previously defined data. Such a model has been already proposed in [START_REF] Pfister | A proposed meta-model for formalizing systems engineering knowledge, based on functional architectural patterns[END_REF] to support the systems engineering (SE) processes defined by the International Council on Systems Engineering (INCOSE). Unfortunately, this proposal considers only the normal, faultless, operation of the system and cannot be used as it is for critical systems.

Several worthwhile recent results may be considered to remove this limitation, however. A method to integrate two classical fault forecasting methods: Failure Mode Effects and Criticality Analysis (FMECA) and Fault Tree Analysis (FTA) in SE processes, is proposed in [START_REF] Guillerm | Combining FMECA and Fault Trees for declining safety requirements of complex systems[END_REF]. An UML profile, termed SOPHIA, for integrating risk analysis in these processes is described in [START_REF] Cancila | SOPHIA: a modeling language for model-based safety engineering[END_REF] while [START_REF] Michael | Hazard analysis and validation metrics framework for system of systems software safety[END_REF] presents a framework for hazard analysis of systems of systems software. Last, the dysfunctional behavior database defined in [START_REF] David | Reliability study of complex physical systems using SysML[END_REF] allows the dysfunctional behaviors be considered through a relevant refinement of the failure mode concept, for physical systems with non-repairable components.

Nevertheless, these valuable results are not fully appropriate to deal with the class of systems which are considered in this work, which focuses on engineering of power production systems [START_REF] Billinton | Reliability Evaluation of Power Systems[END_REF], for the following reasons:

• First of all, since a power plant is built for several decades, its components must be repairable. This constraint is not commonly taken into account in most of dependability analysis where only non-repairable components are considered. The rare authors who discussed this issue ( [START_REF] Raiteri | Repairable fault tree for the automatic evaluation of repair policies[END_REF] and [START_REF] Chaux | Towards a unified definition of Minimal Cut Sequences[END_REF]) do not integrate their work into a SE process. • Dependability analysis assumes very often that the objective of the system is fixed, which is no more true for power plants (and for numerous other critical systems: airplanes, chemical processes, etc.) which are phased mission systems [START_REF] Burdick | Phased mission analysis: A review of new developments and an application[END_REF].

I

and not merely binary components with one faultless state and one faulty state [START_REF] Levitin | Redundancy Optimization for Series-Parallel Multi-State Systems[END_REF]. • This last feature implies that several redundancy policies are possible according to the current mission phase. A faulty component may be replaced by another component which was previously inactive in a given phase and by an already active component which changes its operation mode (speed increase for instance) in another phase. This non-usual kind of redundancy policy may be seen as a resilience strategy [START_REF] Haimes | On the Definition of Resilience in Systems[END_REF]. As depicted on Figure 1, the aim of this paper is to propose a meta-model that contains all the relevant concepts for the class of systems described above and that shall be connected to a meta-model to support SE processes. The result will permit to support processes where functional and dysfunctional analysis will cooperate in a seamless manner to assess dependability attributes, like reliability or unavailability.

The construction of this meta-model is addressed in the next section. The merit of this proposal is illustrated by instantiating, in the third section, the meta-model for a small phased mission system: a part of the water supply system of the steam generator of a power plant, and by building from this instantiated model, in the fourth section, a dysfunctional model, in the form of Markov Chains, for unavailability assessment. Finally, concluding remarks and some outlooks are drawn up.

II. META-MODEL CONSTRUCTION

To address the problem of dependability analysis integration into Systems Engineering processes, this article proposes to extend the SE knowledge meta-model defined in [START_REF] Pfister | A proposed meta-model for formalizing systems engineering knowledge, based on functional architectural patterns[END_REF]. This meta-model has been designed to be an aid for building models that comply with the SE processes suggested by the International Council on Systems Engineering (INCOSE). It includes several classes (Context, Need, Requirement, FunctionalArchitecture, PhysicalArchitecture, Interface, etc.) and relations that can be instantiated to describe the main features of a specific system with faultless components. Due to space limitations, it is not possible to show completely this meta-model but a part of it is depicted by Figure 2.

It expresses that a system is composed of an organized set of components and performs an organized set of functions. Functions are performed by components to which they are allocated.

The meta-model proposed in this paper is aiming at extending this result by adding the semantics required to perform dependability analysis on phased mission systems with repairable multi-state components. It is completed by a list of modeling constraints and definitions. The meta-model is represented by using UML class diagrams [15] and the modeling constraints and definitions are expressed in natural language and in OCL (Object Constraint Language [16]).

A. Modeling Phased Mission Systems

As shown at Figure 3, a phased mission system is characterized by several phases. The system structure, failure and recovery modes, or success criteria can change from one phase to another one ( [START_REF] Burdick | Phased mission analysis: A review of new developments and an application[END_REF] and [START_REF] Meshkat | An overview of the phase-modular fault tree approach to phased mission system analysis[END_REF]). Components and functions are not used similarly during the different phases. Characterizing a phase consists in instantiating the links Phase-Component and Phase-Function. Indeed, these links permit to specify respectively which components must be used and which functions must be performed for each phase.

To simplify the representation, some attributes and links of the initial meta-model (Figure 2) are removed in Figure 3; they must be considered nevertheless in the global model. The current phase can be updated by the method setPhase. Therefore, a component state is a pair built with one operation mode and one failure mode. As depicted on Figure 4, the possible states of a component are defined by instantiating its failure and operation modes.

An initial state of the component must be also defined as stated by the Constraint 2 which guarantees its uniqueness. The link between the classes Phase and Operation mode specify the nominal operation mode of a component for the considered mission phase. The existence and uniqueness of the nominal operation mode of each component involved in a given phase is ensured by the Constraint 3.

C. Modeling the Effects of Component States on Function Achievement

This subsection introduces the new class Effect and enhances the definition of the class Function by adding new attributes and methods (Figure 5). The aim of the new class is to relate the states of the components that are allocated to a function to the complete or partial achievement of this function. To meet this objective, a new attribute allocation must be added to the list of attributes of the class Function. The value of this attribute is the set of components which are allocated to the function.

The class Effect models the contribution of a component in a particular state to the achievement of a function to which it is allocated. Since this contribution depends on the active state, an instance of the class Effect must be defined for every As depicted on Figure 5, the attributes of the class Effect are achievementRate and isUnacceptable. The first one quantifies the contribution of the component, in the considered state, to the achievement of the function; in a system where two identical pumps operate in parallel to fill in a tank the achievementRate for a faultless pump is equal to 50% for a faultless pump and 0% for a failed pump, for instance. The second attribute points out the states which correspond to unsafe conditions; the values of this attribute for the different states are defined when instantiating the meta-model but is always False for the states which satisfy Constraint 5. The new attribute goal of the class Function is a threshold on the sum of the achievement rates of the components allocated to a function, i.e. a function will be declared correctly achieved if and only if the current states of the components which are allocated to it provide an overall achievement rate greater than this threshold. This attribute may be updated when the mission phase changes by the method setGoal given below.

D. Introducing Redundancy Policies

Using redundant components to continue to perform a function despite of the failure of other components that were allocated to this function is a well-known and widespread strategy to increase dependability [START_REF] Villemeur | Reliability, Availability, Maintainability and Safety Assessment, Methods and Techniques[END_REF]. This solution implies that the set of components allocated to a function changes during operation and that redundancy policies which specify these changes are defined. Therefore, a new class Redundancy policy as well as new relations must be introduced in the metamodel and two methods added to the class Function (Figure 6). These methods dynamicAllocation and dynamicDeallocation update the allocation attribute of the class Function and are defined below. The new class Redundancy policy and the new relations where this class is involved permit to specify how the faulty components are replaced by redundant components as follows.

Statement 2: If the aimed function is no more satisfied during a given phase, the failed components are deallocated whereas the redundant components whose current state is declared available (neither failed nor already used for another function), are allocated to the aimed function and their current state is set to spare.

It must be noted that, when several redundancy policies are possible after a component failure, priorities over these policies must be defined. Furthermore, the following two constraints are to be satisfied to ensure consistency of the model. Figure 7 shows the complete meta-model built step by step in this section. This model fits the domain needs since it is possible to construct, from this class diagram, object diagrams for phased mission critical systems with repairable multi-state components and several redundancy policies, as it will be illustrated in the next section.

III. BUILDING OBJECT DIAGRAMS FROM THE META-MODEL

The example considered in this section comes from a power plant: it is a part of the water level control system of the steam generator (Figure 8). This system has been previously described in [START_REF] Kothare | Level control in the steam generator of a nuclear power plant[END_REF] and [START_REF] Zhang | Dynamic reliability: Towards efficient simulation of the availability of a feedwater control system[END_REF], in particular to illustrate contributions in dynamic reliability assessment in the latter reference.

Within this system, only the sub-system composed of the two feeding turbo pumps FTP1 and FTP2, darkened part of Figure 8, will be considered in what follows, for space reasons. This sub-system has to perform only one function F: To supply enough water to the steam generator. The pumps FTP1 and FTP2 may fail and be repaired. To increase dependability of the function, the operation mode of each pump must be managed dynamically according to redundancy policies which will be described later.

Three instance diagrams completed by tables will be necessary to model this system. It must be underlined that the instance diagrams for a larger system must be constructed with the aid of a dedicated software tool, like for instance arKItect © [START_REF] Aboutaleb | An integrated approach to implement system engineering and safety engineering processes: Sasha project[END_REF].

A. Defining the Mission Phases

The main mission of the plant is to produce electric power and is decomposed in three phases (Table I). Function F is mandatory in all phases but only one pump is necessary to perform this function during first and third phases.

B. Defining the Component States

According to the first instance diagram (Figure 9), each pump has three operation modes: OFF (as every component), Run and Overspeed 1 . It will be assumed that FTP1 is the main pump and FTP2 a spare pump for the phases P1 and P3; hence, the initial modes of FTP1 and FTP2 are respectively Run and OFF. The pumps are considered faultless at the initial state; therefore the initial failure mode is OK for both and each pump may fail according two failure modes called Leak and Rupture.

Table II contains the values of the attributes (failureRate, repairRate) of the State class instances 2 . In this study, the probability rates are assumed to be constant for a couple (operation mode, failure mode). The first value (e.g. 0.01 for the state Run-Leak) is the failure rate of the transition that leads to this state and the second one (0.1 for the same state), the repair rate to leave this state. These attributes are not relevant for the states of the first column which are faultless. It must be noted that a pump can fail only if it is active (with a higher probability in the 1 Overspeed means that the speed of the pump is not the nominal value but a higher speed that fits the physical limitations of this component however. 2 These values are rounded average values obtained from several tens of power plants during about forty years. operation mode Overspeed than in Run) 3 . Moreover, the Leak and the Rupture can be repaired if the pump is disabled whereas only the Leak can be repaired if the pump is in the operation mode Run (even if the repair time is longer than in OFF), and no repair is possible if the pump is in the operation mode Overspeed. As the pumps are assumed identical, only one table is necessary.

C. Describing the Effects of the Component States

The contribution of one pump to the function depends on its state. When the pump is disabled (operation mode OFF), this contribution is obviously equal to zero. This is also the case when the Rupture failure has occurred. The analysis is not so simple for the remaining four states: Run-OK, Run-Leak, Overspeed-Ok and Overspeed-Leak.

The attribute goal of function F must be first defined. For phased mission systems this goal is usual equal to 100% for the most demanding phase, P2 in this study, and to smaller values for the other phases. Expert knowledge is mandatory to set these values; it will be assumed hereafter that the goal of F during the phases P1 and P3 is equal to 60%, i.e. that 60% of the maximal water flow is sufficient during these two phases.

Once the goal set for each phase, it is possible to define the values of the attribute achievementRate for every state of a pump (Table III). According to this table, when the operation mode is Run, function F is achieved at 60% when the pump is faultless and 50% when it is leaking. Hence, only one faultless pump is necessary to meet the goal of F during the phase P1 and P3 whereas two running pumps, faultless or leaking, are mandatory for the second phase. Last, a faultless over-speeded pump is sufficient to meet the goal, whatever the phase. This table will permit to define the redundancy policies in the next section. Moreover it is assumed that the states built over the failure mode Rupture are forbidden by taking into account safety. Then, the attribute isUnacceptable of each instance of the class Effect linked to these states is True.

D. Describing the effects of the component states

Three redundancy policies are defined to improve the dependability of the system. Hence, the class Redundancy policy is instantiated three times:

• R1a is defined for the phases P1 and P3, and consists in replacing the pump FTP1 by FTP2, if function F is not satisfied and FTP2 is available. • R1b is defined for the phases P1 and P3, and consists in replacing the pump FTP2 by FTP1, if function F is not satisfied and FTP2 is available. • R2 is defined for the phase P2, and consists in forcing the operation mode of the faultless pump to Overspeed if function F is not satisfied because the other pump has failed.

It can be noted that, during the second phase, if the two pumps fail with the failure mode Leak, according to the table III, the function is correctly achieved anyway and this failure mode can be repaired in the operation mode Run. Then the Leak failure mode does not trigger the redundancy policy R2.

The instance diagram that represents graphically these policies and their relations to other class instances has been split in Figures 10 and11, for clarity reasons.

The attribute priority is not specified because there is no concurrency between these redundancy policies. It can be noted that for the redundancy policy R2 the two pumps are both defined failed and redundant, because they could play the two roles: either FTP1 has failed and FTP2 is redundant, or FTP2 has failed and FTP1 is redundant.

The system features during each phase are summarized in table IV. IV. USING THE META-MODEL FOR DEPENDABILITY ANALYSIS It has been claimed in the introduction that one benefit of a meta-model is to ease the automatic construction of some specific models from previously defined data. This claim will be illustrated in this section where a method to construct systematically the continuous Markov Chains (MC) that represent the behavior of a phased mission system which includes redundant components will be proposed. Continuous MCs are indeed common and relevant models to assess dependability attributes, like (un)availability; systematic construction of such models from instance diagrams is surely a positive consequence of the integration of dependability concerns into a meta-model for SE processes. This construction method is described and exemplified, on the basis of the previous example, in the first sub-section, while analytic assessment of unavailability, from these chains, is dealt with in the second sub-section. It must be underlined that, when dealing with larger systems, both steps must be automated. Automatic construction of the MC will be based on the systematic procedure described below. The reader interested in automatic assessment of unavailability from large MC models will find more details on this issue in [START_REF] Stewart | Introduction to the Numerical Solution of Markov Chains[END_REF].

A. Systematic Construction of MCs from Instance Diagrams

Modeling the evolutions of a phased mission system with repairable components by MCs is a classical approach in the field of dependability analysis ( [START_REF] Alam | Quantitative reliability evaluation of repairable phased-mission systems using Markov approach[END_REF], [START_REF] Murphy | Reliability analysis of phased-mission system: A correct approach[END_REF]). One MC must be constructed for each phase and a transition matrix specifies how the probabilities are distributed over the states when the active phase changes.

From the instance diagrams presented in section III, it is possible to build systematically the MCs depicted in Figures 12 and13; the MC of Figure 12 describes the system behavior for both phases 1 and 3 whereas this behavior for the second phase is given at Figure 13. These MCs have been built systematically, to avoid errors that may appear happen when no method is employed, according to the following generic procedure:

1) The initial state of the chain is defined from the initial states of the components of the system, for the considered phase. 2) For each state of the chain, an outgoing transition is added for each possible failure or repair event. 3) If this transition does not lead to a state which has been already defined, a new state is introduced by taking into account the redundancy policies. 4) The system is unavailable if the sum of the achievement rates of the components in this state is smaller than the goal of the function. This procedure stops when it is no longer possible to add any transition or state.

The MC of Figure 12 has been built by using the abovedefined procedure:

1) The initial states of the pumps are respectively Run-OK and OFF-OK for the phases 1 and 3 (table IV).

Then, the initial state of the chain gathers these two states. 2) Two transitions are starting from the initial state of the chain to represent the two possible failure modes Leak and Rupture. The failure rates are given by the first term of the cells Run-Leak and Run-Rupt in the table II.

3) The states 2 and 3 are introduced to represent a system where one failure has occurred in one component and the other component has been activated according to the appropriate redundancy policy. As the failure mode Rupture is unacceptable, the component for which this failure occurred is forced to the operation mode OFF (state 3). 2) From these new states, the failed pump can be repaired with repair rates given by the second term of the cells Run-Leak and OFF-Rupt in Table II, what leads to the initial state, or the running pump can fail with failure rates given by the first term of the cells Run-Leak and Run-Rupt in the table II.

3) The states 4, 5 and 6 are then introduced to model a system where failures have occurred in the two components. 2) From these states, the failed pumps can be repaired, what leads to the already known states 2 or 3. It shall be noted that the probability rate associated to the transition from the state 4 to the state 2 (respectively 6 to 3) is two times the repair rate defined in the cell Run-Leak (respectively OFF-Rupt) because the two pumps are in the same state; therefore the probability to have a pump repaired at a given time, assuming that they can be repaired simultaneously, is twice better. 4) As the goal of the function in these phases is equal to 60.0, the states 5 and 6 correspond to an unavailable system because the sum of the achievement rates of the two pumps is indeed equal to 50.0 for the state 5 and 0.0 for the state 6, according to the table III.

The MC for the second phase (Figure 13) is built in the same fashion. Two faultless pumps are running in the state 1 and the transition from this state to the state 3 models the redundancy policy R2. The states 2 and 4 are identical to those of Figure 12 but it shall be underlined that no redundancy policy is involved to define these states in the second phase; the transitions between the states 2 and 4 are merely provoked by failure or repair events. The state 5 differs from the corresponding state in Figure 12 because the active pump is over-speeded; there is no transition from this state to the state 3 because the Leak failure cannot be repaired when the operation mode is Overspeed. Last, some transition rates are different because either two pumps are in the same operation mode in the source state of the transition (transition from 1 to 2, 1 to 3, 4 to 2, 6 to 3) or because the rate defined in the table II is different. The darkened states (5 and 6) correspond to an unavailable system; the sum of the achievement rates is smaller than the goal.

The model of the whole behavior of the system for the three phases is given at Figure 14. This model is composed of the three previous MCs as well as three transition matrices φ i→j that allow the computation of the probability to be in a state l of the MC for the phase j at the date t ij , where t ij is the date of a phase change from the phase i to the phase j, from the probabilities to be in a state k of the MC for the phase i at the date t ij . In our case, these matrices are identity matrices, i.e.

∀ ∈ ⟦0,6⟧ � 1 ( 12 ) = 2 ( 12 ) 2 ( 23 ) = 3 ( 23 ) 3 ( 31 ) = 1 ( 31 ) (1)

B. Unavailability Assessment

Once this model obtained, it is possible to compute the unavailability of the system by analytic calculus or Monte-Carlo simulation. The first technique was selected for this work because it can be easily applied to small-sized systems and provides generic solutions (more details on this issue may be found in [START_REF] Stewart | Introduction to the Numerical Solution of Markov Chains[END_REF]).

For the phases 1 and 3, the transition matrix of the MC is (cf. Figure 12):

13 = ⎣ ⎢ ⎢ ⎢ ⎢ ⎡ -0.011 0.1 0.1 0 0 0 0.01 -0.111 0 0.2 0.1 0 0.001 0 -0.111 0 0.1 0.2 0 0.01 0 -0.2 0 0 0 0.001 0.01 0 -0.2 0 0 0 0.001 0 0 -0.2 ⎦ ⎥ ⎥ ⎥ ⎥ ⎤ Let 13 () = [ 1 13 (), 2 13
(), 3 13 (), 4 13 (), 5 13 (), 6 13 ()],

be the probability to be in each state of the chain at the date , 13 is solution of the system of differential equations:

13 () = 13 (). 13

It can be expressed as:

13 () = 13 (0). 13 From this result, the asymptotic unavailability for each one of the phases 1 and 3 can be computed from the asymptotic probabilities to be in the states 5 and 6 of the chain: The asymptotic unavailability for the phase 2 can be calculated in a similar manner: ̅ 2 (∞) = 7.156 × 10 -3 .

The result obtained for a sequence of twelve identical missions where the first phase lasts one day, the second 28 days and the third also 1 day is given at Figure 15. The average unavailability is equal to 0.623827%, which corresponds to approximately 2 days and 6 hours in absolute value.

V. CONCLUSIONS AND OUTLOOKS

This paper has proposed a meta-model based on [START_REF] Pfister | A proposed meta-model for formalizing systems engineering knowledge, based on functional architectural patterns[END_REF] to support the integration of dependability analysis into system engineering processes. This meta-model has been defined for a wide class of systems: phased mission systems with multistate and repairable components and several redundancy policies; it can be easily specialized to deal with more restricted classes of systems, like, for instance, systems with binary and non-repairable components or systems whose mission comprises only one phase.

The benefit of this proposal has been shown on a small example of a critical system coming from a power plant. A systematic procedure to construct MCs which are classical models for availability analysis, from object diagrams derived from the meta-model, has been defined and exemplified.

On-going works are aiming to automate the construction of object diagrams and MCs to be able to deal with full-scale industrial cases. Automatic construction of MCs will take into account uncertainties on the failure and repair rates to better assess the impact of these uncertainties on the final results (probabilities to be in a given state, (un)availability).

Introduction of new classes, attributes and relations in the meta-model, to integrate the new concepts recently presented in the domain of dynamic dependability analysis ( [START_REF] Chaux | Towards a unified definition of Minimal Cut Sequences[END_REF], [START_REF] Rauzy | Sequence algebra, sequence decision diagrams and dynamic fault trees[END_REF]), is also under investigation. This will allow, in particular, to integrate faulty behaviors due to software or hardware components of the Instrumentation and Control system, into the dependability analysis. Finally, validation of this meta-model on case studies larger than the small example which has been presented in this paper will permit to really assess its usefulness, completeness, and flexibility and to estimate correctness and scalability of the analysis methods of the derived models.

Fig. 1 .Fig. 2 .Fig. 3 .

 123 Fig. 1. Contribution of the article

Constraint 2 :

 2 Every component must have one unique initial state. context Component inv: self.operation mode ->one(om: Operation mode | om.isInit = True) and self.failure mode ->one(fm: Failure mode | fm.isInit = True)

Constraint 3 :

 3 A mission phase must define the nominal operation mode for each component which is involved in this phase. context Phase inv: self.Component.Operation mode ->one(om: Operation mode | om = self.nominal) Moreover, the class diagram represented at Figure 4 permits to model the stochastic evolutions of the component in the form of transitions from a state (om i , fm i ) to a state (om i , fm j ). These transitions are provoked by failure and repair events and are defined as follows (Statement 1). Statement 1: Let C be a component, OK be its non-faulty failure mode, and (om,fm) be a faulty state of C where the attributes failureRate and repairRate are not null: • If C is in the faultless state (om,OK), then it can fail according to the failure mode fm and the transition from (om,OK) to (om,fm) occurs with the specified failure rate. • If C is in the faulty state (om,fm), then it can be repaired and the transition from (om,fm) to (om,OK) occurs with the specified repair rate. The current state of a component can be updated by the method setState defined below. Definition 2: Method setState. context Component::setState(P: Phase) pre: self.Operation mode.State ->includes(S) post: self.currentState = S

Fig. 4 . 4 :

 44 Fig. 4. Step 2: definition of the component states

Constraint 5 :

 5 A disabled (operation mode OFF) or faultless (failure mode OK) component has no unacceptable effect. context Effect inv: self.state.operation mode.name = 'OFF' or self.state.failure mode.name = 'OK' implies self.isUnacceptable = False

Definition 5 :

 5 Method dynamicAllocation. context Function::dynamicAllocation(C: Component) pre: self.static allocation ->includes(C) post: self.allocation = self.allocation@pre ->including(C) Definition 6: Method dynamicDeallocation. context Function::dynamicDeallocation(C: Component) pre: self.static allocation ->includes(C) post:self.allocation = self.allocation@pre ->excluding(C)

Fig. 5 .

 5 Fig. 5. Step 3: definition of the effects Fig. 6. Step 4: definition of the redundancy policies

Constraint 6 :

 6 The states whose roles are available or spare for a redundancy policy must be linked to the redundant components of this redundancy policy. context Redundancy policy inv: self.redundant.operation mode.state ->includesAll (self.available ->union(self.spare)) Constraint 7: Exactly one spare state must be defined for each redundant component. context Redundancy policy inv: self.redundant ->forAll(c: Component | c.Operation mode.State ->one(s: State | self.spare ->includes(s)))

Fig. 7 .Fig. 8 .

 78 Fig. 7. Complete meta-model to support the integration of dependability analysis into SE processes

Fig. 10 .

 10 Fig. 10. Second partial instance diagram for the considered sub-system

Fig. 12 . 3 Fig. 13 .

 12313 Fig. 12. MC for the phases 1 and 3

̅ 13 (

 13 ∞) = 513 (∞) + 613 (∞) = 9.41 × 10 -4

Fig. 14 .3→1Fig. 15 .

 1415 Fig. 14. MCs for a mission

Definition 1 :

 1 Method setPhase. Every component must have at least one operation mode noted OFF and one failure mode noted OK.

	context System::setPhase(P: Phase) pre: self.Phase ->includes(P)
	post: self.currentPhase = P
	B. Modeling the Component States
	Each component can be activated and fail according to
	several operation and failure modes. These modes represent
	respectively the functional and dysfunctional properties of the
	component. Moreover, at least one operation mode: inactive
	(noted OFF), and one failure mode: faultless (noted OK) must
	exist (Constraint 1).
	Constraint 1:

context Component inv: self

.operation mode ->one (om | om.name = 'OFF') and self.failure mode ->one(fm | fm.name = 'OK')

  self.goal = newValue Last, the method isSatisfied checks whether the sum of the achievement rates of the components allocated to a function, in their current states, is greater than goal.

	Definition 4: Method isSatisfied().
	context Function::isSatisfied():boolean body: self.effect ->select(e:Effect |
	e.state.operation mode.component.currentState = e.state and self.allocation ->includes(
	e.state.operation mode.component)).achievementRate ->sum() ≥ self.goal

Definition 3: Method setGoal. context Function::setGoal(newValue: float) pre: newValue ≥ 0.0 post:

TABLE IV SUMMARY

 IV OF THE FEATURES OF THE INSTANCES OF THE CLASS Phase

	Phase	Nominal OM FTP1 FTP2	Function goal	Redundancy policies
	P1	Run	OFF	60.0	R1a, R1b
	P2	Run	Run	100.0	R2
	P3	Run	OFF	60.0	R1a, R1b

The state OFF-Leak (respectively OFF-Rupture) is then not reachable from the state OFF-OK but from the state Run-Leak (respectively Run-Rupture). Fig. 9. First partial instance diagram for the considered sub-system


Visual cortex as a depth computer: Grossberg's LAMINART or how we see the world in depth

Birgitta Dresp-Langley ICube, UMR 7357 CNRS-UdS, Strasbourg, FRANCE

The laminar structure of cortex was first described by [START_REF] Ramón | Comparative study of the sensory areas of the human cortex[END_REF]. It consists of six different layers, each containing a characteristic distribution of neurons and cell types connecting with other cortical and sub-cortical regions. These extensive connections form microcircuits that group into functionally distinct columns (e.g. [START_REF] Mountcastle | The columnar organization of the neocortex[END_REF]. To account for the Venetian Blind Effect, Cao and Grossberg here propose a model which exploits the functional organization of visual cortex, in particular the selective firing activities of visual cortical neurons responding to the orientation and/or the luminance contrast of surface boundaries, demonstrated first by Hubel &[START_REF] Hubel | Receptive fields of single neurons in the cat's striate cortex[END_REF][START_REF] Hubel | Receptive fields, binocular interaction, and functional architecture in the cat's visual cortex[END_REF]. Orientation selective neurons in visual cortex belong to an array of orientation columns which, combined with the eye of origin, provide a basic representation for visual information processing. The resulting hypercolumns are not functionally independent, but display long-range horizontal functional interactions, generating a seemingly infinite number of possibilities for information processing.

Cao and Grossberg's approach exploits these possibilities by suggesting a computational model architecture where selective cortical signals, activated by boundaries and the luminance contrasts which fill in the surfaces they delineate, control what the authors call three-dimensional surface capture. Rather than looking for some unique binocular matching rule to account for the effect considered here (the Venetian blind effect), they suggest a cortical disparity filter, which combines both monocular and binocular information to generate visual depth. On the basis of combined monocular and binocular input, the filter selectively activates the most pertinent perceptual groupings for generating the most likely depth percepts. It suppresses other, less likely but possible, groupings, and the winner takes all in this race for structural meaning, where cooperative and competitive interactions determine the relative synaptic weights of the different ongoing processes.

Three-dimensional surface depth is thus predicted to be a result of complex, facilitating and suppressing, layer-to-layer neural interactions. The mechanistic dynamics of this model have the functional advantage of being closely linked to physiological evidence for long-range feed-back interactions almost everywhere in cortex, and in particular in the visual cortex (see the references provided by the authors). Such interactions imply that any explanation of a complex perceptual phenomenon in terms of a single process or mechanism is unlikely. As a consequence, the model by Cao and Grossberg works out how monocular input is taken into account in the computation of visual depth at the cortical levels and delivers a complex yet highly plausible functional assumption for monocular depth effects. Monocular cues to pictorial depth, sometimes called monocular depth effects, have been listed by Leonardo da Vinci in his Trattato della Pittura well before they were measured psychophysically by visual scientists. Luminance contrast, relative size, and relative height in the visual field are three such cues, and any of them is self-sufficient in determining what humans may perceive as nearer or further away in a visual image [START_REF] Dresp | On illusory contours and their functional significance[END_REF][START_REF] O'shea | Contrast as a depth cue[END_REF][START_REF] Guibal | Interaction of color and geometric cues in depth perception: When does "red" mean "near[END_REF][START_REF] Dresp-Langley | Color and figure-ground: From signals to qualia[END_REF], for a review).

Monocular processing is essential for understanding the mechanisms that allow oneeyed individuals to see depth. Results of recent psychophysical research lead to suggest that the role of monocular vision in depth perception is, indeed, far more important than previously considered [START_REF] Vishwanathan | Seeing in 3D with just one eye[END_REF]. Finally, by exploiting the laminar functional organization of visual cortex to account for how we see the Venetian Blind Effect (and other phenomena) reconnects visual science with an important aspect of brain evolution. Laminar functional organization of cortex is found in various species and it is likely to have evolved as a result of pressure on the brain's functional anatomy to adapt to constantly changing, visual and other, environments. Seeing in depth is even more relevant to animals as than it is to humans. Animal color vision in particular, which is equipped with anatomical and functional properties beyond those of human perception, plays an important part in how the different species perceive three-dimensional biological reality (e.g. [START_REF] Nolan | Immunocompetence is signalled by ornamental color in King Penguins, Aptenodytes Patagonicus[END_REF]. Many of the phenomena accounted for by Cao and Grossberg's model have been discussed previously in terms of hypotheses relative to the adaptive function of visual processing, as for example the illusory contours seen in the Kanzisa square (e.g. see [START_REF] Dresp | On illusory contours and their functional significance[END_REF], for a review), which the authors here repeatedly use to illustrate some of their model assumptions for the Venetian Blind Effect.


Gergely Acs
email: gergely.acs@inria.fr

Claude Castelluccia
email: claude.castelluccia@inria.fr

A Case Study: Privacy Preserving Release of Spatio-temporal Density in Paris

With billions of handsets in use worldwide, the quantity of mobility data is gigantic. When aggregated they can help understand complex processes, such as the spread viruses, and built better transportation systems, prevent traffic congestion. While the benefits provided by these datasets are indisputable, they unfortunately pose a considerable threat to location privacy.

In this paper, we present a new anonymization scheme to release the spatio-temporal density of Paris, in France, i.e., the number of individuals in 989 different areas of the city released every hour over a whole week. The density is computed from a call-data-record (CDR) dataset, provided by the French Telecom operator Orange, containing the CDR of roughly 2 million users over one week. Our scheme is differential private, and hence, provides provable privacy guarantee to each individual in the dataset. Our main goal with this case study is to show that, even with large dimensional sensitive data, differential privacy can provide practical utility with meaningful privacy guarantee, if the anonymization scheme is carefully designed. This work is part of the national project XData (http://xdata.fr) that aims at combining large (anonymized) datasets provided by different service providers (telecom, electricity, water management, postal service, etc.).

INTRODUCTION

Mobile phone datasets have become widely available in recent years and have opened the possibility to improve our understanding of large-scale social networks by investigating how people exchange information, interact, and develop social interactions. While the benefits provided by these datasets are indisputable, they unfortunately pose a considerable threat to location privacy. Not only this can impact people lives negatively, this also affects research. Because privacy is so important to people, companies and researchers are reluctant to publish mobile phone datasets by fear of being held responsible for potential privacy breaches. It is therefore urgent to develop practical tools for private releases and analysis of mobility datasets.

This paper focuses on applications that only need location counts, i.e., the repartition of visitors on a map at a given period of time. This information can typically be published by dividing a map into cells, and then releasing the count (i.e., number of users) associated to each of the cells. However, it is important to ensure that this publication does not leak any information about the mobility patterns of individual users, which might be trivial when, for example, the count values are small.

Why Differential Privacy? Privacy has different definitions and different models have been proposed. In this paper, we use differential privacy [START_REF] Dwork | Calibrating noise to sensitivity in private data analysis[END_REF] which has emerged as a compelling privacy model. The advantage of this model, compared to the many others proposed in the literature, is two-fold. First, it provides a formal and measurable privacy guarantee regardless what other background information or sophisticated inference technique the adversary uses even in the future. Second, it is closed with respect to sequential and parallel composition, i.e., the result of the sequential or parallel combination of two differential private algorithms is also differential private.

This has particular importance in practice, since it does not only simplify the design of anonymization solutions, but also allows to measure how much privacy remains when a given dataset, or a set of correlated datasets, is anonymized (and released) several times, possibly by different entities.

Differential-private schemes often requires adding noise to the published data (e.g., to the published location counts) so that it leaks almost no information about any participating individual, but still reveals aggregated information about the population as a whole. The variance of the noise is calibrated to both the sensitivity of the counts (i.e., the maximal change of the counts due to the inclusion/exclusion of a single record in a dataset) and a desired privacy level ε. For large-dimensional data, such as temporal density, the sensitivity is usually so large that the added noise is much larger than the actual density count values for stringent privacy requirement (i.e. ε < 1). Consequently, the resultant anonymized data is often meaningless.

Contributions. In this paper, we show that, for a given privacy level (i.e., a given ε), the magnitude of noise can be substantially reduced by using several optimizations and by customizing the anonymization mechanisms to the public characteristics of datasets and applications. We observe that the temporal density, of each cell, can be characterized by a periodic time series. This is explained by the fact that aggregated mobility patterns are quite periodic. Moreover these time series follow very similar trends, as most people in nearby cells have similar calling patterns. As a consequence, time series can be compressed by sampling, clustering and low-pass filtering. Sampling and clustering reduce data sensitivity, which results in lower added noise and better performance. We further attenuate the distortion resulting from the compression and perturbation phases via novel optimization and post-processing algorithms. We show experimentally that the achieved performance is quite high and that differential privacy can be practical in real-world applications. However, we believe that there are probably no differential private anonymization techniques that fit all applications, and that anonymization algorithms have to be customized to each application and dataset.

The XData project. This work was done in the context of the French XData project. XData is a national funded project under the "Big Data program", whose goal is to study the benefit of combining and cross-processing different types of datasets provided by various service providers (such as Orange, Electricité de France, La Poste, etc.). However, according to the European Data Protection laws (Directive 95/46/EC), all datasets have to be anonymized, prior crossprocessing, such that data subjects are no longer identifiable. The law does not dictate any specific privacy model, but stipulates that "to determine whether a person is identifiable, account should be taken of all the means likely reasonably to be used either by the controller or by any other person to identify the said person". We believe that the best existing model to achieve this goal is probably differential privacy. This paper shows how to anonymize the mobility data, provided by Orange, under the differential privacy model. In particular, we show how to release density information. Geographical density is useful in many of applications envisioned by the XData projet, such as identifying areas where to install new businesses or build new infrastructures.

RELATED WORK

Several recent studies have demonstrated the privacy risks of releasing location data by re-identifying individuals from geospatial data sets [START_REF] De Montjoye | Unique in the crowd: The privacy bounds of human mobility[END_REF][START_REF] Zang | Anonymization of location data does not work: A large-scale measurement study[END_REF][START_REF] Golle | On the anonymity of home/work location pairs[END_REF]. As a result, a plethora of privacy-preserving techniques have been introduced, however, most of them do not provide any formal privacy guarantee (see [START_REF] Chow | Trajectory privacy in location-based services and data publication[END_REF] and the references therein).

Differential privacy (DP) was first rigorously presented in [START_REF] Dwork | Calibrating noise to sensitivity in private data analysis[END_REF] with the Laplace mechanism (LPA) as a first generic tool to guarantee DP. There exist two relaxations of ε-DP; (ε,δ)-probabilistic DP [START_REF] Machanavajjhala | Privacy: Theory meets practice on the map[END_REF] and (ε, δ)-indistinguishability [START_REF] Dwork | Our data, ourselves: privacy via distributed noise generation[END_REF]. The former guarantees ε-DP with high probability (≥ 1-δ), while the latter relaxes the bound of ε-DP.

Location privacy has also been addressed in the context of DP. [START_REF] Chen | Differentially private sequential data publication via variable-length n-grams[END_REF] applies DP to location, and more generally, sequential data release. However, this scheme does not release time information apart from the sequentiality of locations. Another recent work [START_REF] Andrés | Geo-indistinguishability: Differential privacy for location-based systems[END_REF] formalizes location privacy as protecting the users' location within a radius r with a level of privacy that depends on r. This corresponds to a generalized version of DP. They target LBS applications and add noise directly to users' GPS coordinates. Commuting patterns in U.S. were anonymized in the probabilistic DP model in [START_REF] Machanavajjhala | Privacy: Theory meets practice on the map[END_REF]. This scheme has also been deployed in practice within the project called OnTheMap by the U.S. Census Bureau. Other authors [START_REF] Cormode | Differentially private spatial decompositions[END_REF][START_REF] Li | Differentially private grids for geospatial data[END_REF] apply different spatial decomposition techniques to partition the two dimensional domain into cells, and then obtain noisy counts for each cell. However, these techniques are not concerned with releasing multiple counts over time.

Several DP techniques have been proposed to release time series. In [START_REF] Fan | Real-time aggregate monitoring with differential privacy[END_REF], the authors propose a framework to release real-time aggregate statistics under DP based on filtering and adaptive sampling. Some other proposals [START_REF] Dwork | Differential privacy under continual observation[END_REF][START_REF] Chan | Private and continual release of statistics[END_REF] provide a weaker guarantee on continuous data streams; they provide event-level privacy to protect an event (i.e., one user's presence at a particular time point), rather than the presence of that user. As all these works address the real-time (interactive) release of aggregates, they are usually less accurate than off-line (non-interactive) approaches, which can access the whole time series and build more accurate models for perturbation.

The most related work to ours is [START_REF] Rastogi | Differentially private aggregation of distributed time-series with transformation and encryption[END_REF][START_REF] Kellaris | Practical differential privacy via grouping and smoothing[END_REF][START_REF] Mir | Dp-where: Differentially private modeling of human mobility[END_REF]. All of them address the off-line release of time series under DP. Rastogi and Nath [START_REF] Rastogi | Differentially private aggregation of distributed time-series with transformation and encryption[END_REF] proposed a Discrete Fourier Transform (DFT) based algorithm which guarantees DP by perturbing the discrete Fourier coefficients. This technique was further improved for DP histogram release in [START_REF] Ács | Differentially private histogram publishing through lossy compression[END_REF]. We reuse the improved private DFT algorithm from [START_REF] Ács | Differentially private histogram publishing through lossy compression[END_REF], and further improve its accuracy at the cost of some privacy. In [START_REF] Kellaris | Practical differential privacy via grouping and smoothing[END_REF], the time series are pre-processed by pre-sampling and smoothing them via averaging. These techniques diminish the global sensitivity of the data, and thereby allows to lower the injected Laplace noise. We also use a similar sampling technique to [START_REF] Kellaris | Practical differential privacy via grouping and smoothing[END_REF] in order to compress time series. However, [START_REF] Kellaris | Practical differential privacy via grouping and smoothing[END_REF] is a more general solution and targets even non-periodic time series where averaging may be a better perturbation model than DFT. As aggregated location counts are typically periodic, DFT is a more accurate perturbation model for our application. Finally, DP-WHERE [START_REF] Mir | Dp-where: Differentially private modeling of human mobility[END_REF] adds noise to the set of empirical probability distributions which is derived from the original CDR datasets, and samples from these distributions to generate synthetic CDRs. Instead of releasing the whole CDR dataset, we only aim at releasing the temporal density of IRIS cells which is a more specific problem. This mapping between CDRs and IRIS cells influences utility and is not considered in [START_REF] Mir | Dp-where: Differentially private modeling of human mobility[END_REF].

PRELIMINARIES

Differential Privacy

Intuitively, differential privacy [START_REF] Dwork | Calibrating noise to sensitivity in private data analysis[END_REF] (DP) requires that the outcome of any computation be insensitive to the change of a single record in the dataset. Consequently, for a record owner, it means that any privacy breach will not be due to participating in the database.

Definition 1 (Differential Privacy) A privacy mechanism A gives ε-differential privacy if for any database D1 and D2 differing on at most one record, and for any possible output O ∈ Range(A),

e -ε × P r[A(D2) = O] ≤ P r[A(D1) = O] ≤ e ε × P r[A(D2) = O]
where the probability is taken over the randomness of A.

A relaxation of DP is probabilistic-DP [START_REF] Machanavajjhala | Privacy: Theory meets practice on the map[END_REF], where privacy breaches may occur with very small probability.

Definition 2 (Probabilistic Differential Privacy [START_REF] Machanavajjhala | Privacy: Theory meets practice on the map[END_REF]) A privacy mechanism A gives (ε, δ)-probabilistic differential privacy if for any database D1 and D2 differing on at most one record, and for any possible output O ∈ Range(A), we can partition the output space Ω into Ω1 and Ω2 such that (1) for all O ∈ Ω1,

e -ε × P r[A(D2) = O] ≤ P r[A(D1) = O] ≤ e ε × P r[A(D2) = O]
and (2) for any database D, P r[A(D) ∈ Ω2] ≤ δ where the probability is taken over the randomness of A.

The latter definition guarantees that algorithm A achieves DP with high probability (≥ 1 -δ), and the set Ω2 contains all outputs that are privacy breaches according to Definition 1. The probability of these outputs are bounded by δ. Notice that with δ = 0 probabilistic DP boils down to Definition 1. Although probabilistic DP has weaker privacy guarantee than Definition 1, it provides higher utility in practice.

The definition of differential privacy enjoys the property of sequential composition, which specifies the privacy guarantee in a sequence of computation.

Theorem 1 ( [START_REF] Mcsherry | Privacy integrated queries: an extensible platform for privacy-preserving data analysis[END_REF]) Let Ai each provide εi-differential privacy. A sequence of Ai(D) over the dataset D provides i εi-differential privacy.

Differential Private Mechanisms

Three principal techniques for achieving (probabilistic) DP are Laplace mechanism [START_REF] Dwork | Calibrating noise to sensitivity in private data analysis[END_REF] (LPA), Gaussian mechanism [START_REF] Götz | On User Privacy in Personalized Mobile Services[END_REF] (GPA), and Exponential mechanism [START_REF] Mcsherry | Mechanism design via differential privacy[END_REF]. A fundamental concept of all these techniques is the global sensitivity of a function [START_REF] Dwork | Calibrating noise to sensitivity in private data analysis[END_REF]:

Definition 3 (Global Sensitivity) For any function f : D → R d , the sensitivity of f is ∆f = maxD 1 ,D 2 ||f (D1) - f (D2)||1 for all D1, D2 differing in at most one record.
The global sensitivity is also called as L1-sensitivity due to the L1-norm used in its definition and is denoted by ∆1f . Similarly, the L2-sensitivity ∆2f of a function f , which is used later in this paper, is defined by the L2-norm || • ||2.

Laplace mechanism (LPA).

A standard mechanism to achieve differential privacy is to add Laplace noise to the true output of a function. The noise is generated according to a Laplace distribution with the probability density function (PDF) p(x|λ) =1 2λ e -|x|/λ , where λ is calibrated as follows. Theorem 2 ([10]) For any function f :

D → R d , the mechanism A A(D) = f (D) + L1(λ), . . . , L d (λ)
gives ε-differential privacy, if Li(λ) are i.i.d Laplace variables with scale parameter λ = ∆f /ε.

Gaussian mechanism (GPA). An alternative technique

to achieve probabilistic DP is to add Gaussian noise to the true output of a function. The noise is generated according to the Gaussian distribution with the PDF p

(x|σ) = 1 √ 2πσ e -x 2 /2σ 2 . Theorem 3 ([15]) For any function f : D → R d , the mechanism A A(D) = f (D) + G1(σ), . . . , G d (σ)
gives (ε, δ)-probabilistic differential privacy for any ε ≤ 1 and σ2 ≥ 2(∆2f ) 2 ln(4/δ)/ε 2 , where Gi(σ) are i.i.d Gaussian variables with variance σ 2 .

Assuming identical values of ε, a Gaussian random variable is more concentrated around 0 than a Laplace random variable thereby ensuring better utility for GPA. However, this larger accuracy also entails weaker privacy, since there is a small probability δ that ε-DP will not hold.

Exponential mechanism. The exponential mechanism [START_REF] Mcsherry | Mechanism design via differential privacy[END_REF] captures all DP mechanisms with a measurable output space. In particular, it assigns exponentially greater probabilities of being selected to outputs of higher utility so that the final output would be close to the optimum yet still differential private.

Utility Metrics

The error between the private and original time series is measured by the following metrics. Consider counts X = {X0, . . . , Xn-1}. We denote the original time series by X, the sanitized one by X.

Mean Relative Error (MRE

): MRE(X, X) = (1/n) n-1 i=0 | Xi -X i | max(γ,X i )
, where the sanity bound γ mitigates the effect of very small counts.

Following the convention [START_REF] Xiao | iReduct: Differential privacy with reduced relative errors[END_REF], we adjust γ to 0.1% of n-1 i=0 Xi.

Pearson Correlation (PC): PC(X, X) = n-1 i=0 (X i -i X i /n)( Xi -i Xi /n) n-1 i=0 (X i -i X i /n) 2 n-1 i=0 ( Xi -i Xi /n) 2
. PC measures the linear correlation between the noisy and the original time series (i.e., whether they have similar trends), and it always falls between -1 and 1.

PROBLEM DEFINITION

Our goal is to release the spatio-temporal density of 989 non-overlapping areas in Paris, called IRIS cells. Each cell is defined by INSEE 1 and covers about 2000 inhabitants. The set of all IRIS cells is denoted by L henceforth, and are depicted in Figure 1 based on their contours 2 .

We aim to release the number of all individuals who visited a specific IRIS cell in each hour over a whole week. Since human mobility trajectories exhibit a high degree of temporal and spatial regularity [START_REF] Gonzalez | Understanding individual human mobility patterns[END_REF], one week long period should be sufficient for most practical applications. Therefore, we are interested in the time series

X L = X L 0 , X L 1 , . . . , X L 167
of any IRIS cell L ∈ L, where X L t denotes the number of individuals at L in the (t + 1)th hour of the week, such that any single individual can visit a tower only once in an hour.

We will omit t and L in the sequel, if they are unambiguous in the given context. X L denotes the set of time series of all IRIS cells in the sequel.

Dataset

To compute these time series, we use a CDR (Call Data Record) dataset provided by the French telecom company Orange3 , where T represents the set of cell towers of the operator, and a cell tower T ∈ T is visited by an individual at time t, if the operator has a recorded event at time t at tower T related to the individual. An event can be an incoming/outgoing call or message to/from the individual. This dataset contains the events of N = 1, 992, 846 users at |T| = 1303 towers within the administrative region of Paris (i.e., the union of all IRIS cells) over a single week (10/09/2007 -17/09/2007). Within this interval, the average number of events per user is 13.55 with a standard deviation of 18.33 (assuming that an individual can visit any tower cell only once in an hour) and with a maximum at 732. The set of all events related to an individual constitute his/her record (trajectory) in the dataset. Similarly to IRIS cells, we can create another set of time series X T , where X T t denotes the number of visits of tower T in the (t + 1)th hour of the week.

Computing the time series of IRIS cells

We map the counts in X T to X L as follows. First, we compute the Voronoi tesselation of the towers cells T which is shown in Figure 1. Then, we calculate the count of each IRIS cell in each hour from the counts of its overlapping tower cells; each tower cell contributes with a count which is proportional to the size of the overlapping area. More specifically, if an IRIS cell L overlaps with tower cells {T1, T2, . . . , Tc}, then

X L t = c i=1 X T i t × size(Ti ∩ L) size(Ti) (1)
at time t. The rationale behind this mapping is that users are usually registered at the geographically closest tower at any time. We acknowledge that this mapping algorithm might sometimes be incorrect, since the real association of users and towers depends on several other factors such as signal strength or load-balancing. Nevertheless, without more details of the cellular network beyond the towers' GPS position, we are not aware of any better mapping technique.

PRIVACY PRESERVING RELEASE OF SPATIO-TEMPORAL DENSITY

Our aim is to transform the time series of all IRIS cells X L to a sanitized version XL such that XL satisfies Definition 1. That is, the distribution of XL will be insensitive (up to ε) to all the visits of any single user during the whole week, meanwhile the error between XL and X L is small.

Our sanitization algorithm is sketched in Algorithm 1. First, the input dataset is pre-sampled such that only ℓ visits are retained per user (Line 1). This ensures that the global sensitivity of all the time series (i.e., X L ) is no more than ℓ. Then, the pre-sampled time series of each IRIS cell is computed from that of the tower cells using Voronoi-tesselation 4: Perturb X C into XC (Section 5.3) //see Algorithm 2 5: Apply smoothing on XC (Section 5.4) 6: Compute XL from XC using Formula [START_REF] Ács | Differentially private histogram publishing through lossy compression[END_REF] and Formula (1) (Line 2). After the largest cells, which cover the whole city and also have large counts, from T∪L are identified (Line 3), their time series are perturbed to guarantee privacy (Line 4). In order to mitigate the distortion of the previous steps, we apply smoothing on the perturbed time series as a post-processing step (Line 5). Finally, the time series of all IRIS cells are computed from the post-processed time series in C (Line 6).

Pre-sampling

To perturb the time series of all IRIS cells, we first have to compute their sensitivity, i.e., ∆1(X L ). To this end, we first need to calculate the sensitivity of the time series of all tower cells, i.e., ∆1(X T ). Indeed, Formula (1) does not change the L1-sensitivity of tower counts, and hence, ∆1(X T ) = ∆1(X L ).

∆1(X T ) is given by the maximum total number of (tower) visits of a single user in any input dataset. This upper bound must universally hold for all possible input datasets, and is usually on the order of few hundreds; recall that the maximum number of visits per user is 732 in our dataset. This would require excessive noise to be added in the perturbation phase. We instead follow a different approach and divide the whole sanitization process into two main steps. We first perturb a pre-sample of our dataset which better withstands perturbation, and then mitigate the distortion effect of sampling in a post-processing step described later in Section 5. [START_REF] Chan | Private and continual release of statistics[END_REF].

In particular, we truncate each record of any input dataset by considering at most one visit per hour for each user, and then select at most ℓ of such visits per user uniformly at random over the whole week. This implies that a user can contribute with at most ℓ to all the counts in total regardless of the input dataset, and hence, the L1-sensitivity of the dataset always becomes ℓ. The pre-sampled dataset is denoted by X, and ∆1(X T ) = ∆1(X L ) = ℓ.

Computing the largest covering cells

To sanitize X L , there are two basic (naive) approaches.

First, we can directly perturb the IRIS counts X L by applying the Laplace mechanism: XL t = X L t + L(ℓ/ε) for all L ∈ L. Alternatively, we can first perturb the tower counts X

T to obtain XT , then compute the noisy IRIS counts XL from XT by applying Formula (1). Although both techniques guarantee ε-DP according to Theorem 2, in terms of utility, they are suboptimal. Cells with small counts have larger relative error, whereas larger counts better resist noise. This is due to the fact that the injected noise is independent of the magnitude of the original counts but only depends on their sensitivity. Therefore, the best approach is to first select cells having the largest counts (which can be either tower or IRIS cells) such that they cover whole Paris, perturb the counts of these cells, and then recompute the noisy counts of the smaller IRIS cells, which were not selected in the cover, from the larger (noisy) tower counts.

However, selecting the optimal cover of Paris (i.e., the set of cells having the largest counts) must also be differential private. Fortunately, a simple heuristic helps us to accurately approximate the optimal cover without using the true counts of the cells (that would require to introduce more noise): cells with large size tend to have large counts 4 . This is also confirmed by Figure 2a and2b. Hence, we re-state the problem as follows.

How can we select the minimum cardinality subset of cells C ⊆ T ∪ L such that C is a complete cover (i.e., covers whole Paris)?

More formally, let G(V, E) denote a graph, where each vertex corresponds to a cell in T ∪ L, and (v, v ′ ) ∈ E iff cells v and v ′ overlap. In this setting, our problem translates to the classical minimum vertex cover problem [START_REF] Bondy | Graph Theory with Applications[END_REF]. Indeed, as T and L are also complete covers of Paris, each vertex in V has at least one edge (a tower cell is always overlapped by at least one IRIS cell, and vice-verse), and we want to compute the minimum cardinality subset of cells which cover all the overlapping areas between cells. Although the minimum vertex cover problem is NP-hard in general, our covering problem belongs to the special cases which can be efficiently solved.

Theorem 5 (Minimum cardinality cover) Selecting the minimum cardinality subset of cells

C ⊆ T ∪ L such that C is a complete cover can be solved in O(|T||L| |T| + |L|).
Proof. G(V, E) is bipartite, since the IRIS cells as well as the tower cells are partitionings of Paris, i.e., there are no overlapping cells in any of the two sets. Hence, for all (v, v ′ ) ∈ E, either v ∈ T and v ′ ∈ L, or v ′ ∈ T and v ∈ L. For bipartite graphs, the minimum vertex cover problem is equivalent to the maximum matching problem based on König's theorem [START_REF] Bondy | Graph Theory with Applications[END_REF], which can be solved in polynomial time, e.g., with the Hopcroft-Karp algorithm [START_REF] Hopcroft | An n 5/2 algorithm for maximum matchings in bipartite graphs[END_REF] Figure 2c shows the largest IRIS and tower cells covering Paris. We computed the mean count of each cell over the whole week which are illustrated by the cell colors. Apparently, the minimum cover contains cells with larger counts; the mean counts are 91 on average for the IRIS cells (Figure 2b) and 120 on average for the minimum cover (Figure 2c).

However, care must be taken before computing the cell counts in the minimum cover C. Since C can contain both IRIS and tower cells which may overlap, the L1-sensitivity of all the counts in C can be larger than ℓ. Indeed, if C contains a tower cell T and one of its overlapping IRIS cell L, then adding/removing a user who visited T at time t will change X T t with 1, and also X L t with a non-zero value. A trivial (but not optimal) solution is to modify the counts of all towers in C which have overlapping IRIS cells. For example, if T overlaps with IRIS cells {L1, L2, . . . , Lc}, then X T t should be reduced by c i=1 X

L i t × size(L i ∩T ) size(T ) .

Perturbation

After identifying the largest covering cells, their time series (i.e., X C ) can be perturbed by adding L(ℓ/ε) to each count in all time series (see Theorem 2). Unfortunately, this naive method provides very poor results which is also illustrated in Figure 3a. Indeed, individual cells have much smaller counts than the magnitude of the injected noise; the standard deviation of the Laplacian noise is 141 with ε = 0.3, which is even larger than the mean count in the minimum cover.

A better approach exploits (1) the similarity of geographically close time series, as well as (2) their periodic nature. In particular, we first cluster nearby less populated cells until their aggregated counts become sufficiently large to resist noise. The key observation is that the time series of close cells follow very similar trends, but their counts usually have different magnitudes. Hence, if we simply aggregate (i.e., sum up) all time series within such a cluster, the aggregated series will have a trend close to its individual components yet large enough counts to tolerate perturbation. To this end, we first accurately approximate the time series of individual cells by normalizing their aggregated time series (i.e., divide the aggregated count of each hour with the total number of visits inside the cluster), and scale back with the (noisy) total number of visits of individual cells.

In order to guarantee DP, we also need to perturb the aggregated time series before normalization. To do so, we exploit their periodic nature and apply a Fourier-based perturbation scheme [START_REF] Rastogi | Differentially private aggregation of distributed time-series with transformation and encryption[END_REF][START_REF] Ács | Differentially private histogram publishing through lossy compression[END_REF]: we add noise to the Fourier coefficients of the aggregated time series, and remove all high-frequency components that would be suppressed by the noise. As only low-frequency components are retained and perturbed, this method preserves the trends of the original data more faithfully than LPA.

The whole perturbation process is summarized in Algorithm 2. First, the noisy total number of visits of each cell in the minimum cover C is computed by adding noise L(2ℓ/ε) to 167 t=0 X i t for cell i (Line 1). These noisy total counts are used to cluster similar cells in Line 2 by invoking Algorithm 3. When the clusters are created, their aggregated time series (i.e., the sum of all cells' time series within the cluster) is perturbed with a Fourier-based perturbation scheme in Line 5 (Algorithm 4). Finally, the perturbed time series of each cell i in cover C is computed in Line 7 by scaling back 2: E := Cluster(C, τ, Ŝ) //see Algorithm 3 3:

for each cluster E ∈ E do 4: X E := i∈E X i 0 , i∈E X i 1 , . . . , i∈E X i 167

5:

XE := EFPAG(X E , ε/2, δ) //see Algorithm 4 6: for each cell i ∈ E do 7:

Xi := Ŝi • ( XE t /|| XE || 1 )
8: end for 9: end for the normalized aggregated time series with the noisy total count cell i (i.e., with Ŝi). Since Line 1 guarantees ε/2-DP to the total counts (∆1(X C ) = ℓ), it follows from Theorem 1 that Algorithm 2 is (ε, δ)-DP if EFPAG is (ε/2, δ)-DP.

Clustering cells

Algorithm 3 is a simple iterative process that, in each iteration, merges the least visited cluster with its geographically closest neighboring cluster until all clusters in the resultant configuration have a total count larger than a predefined threshold τ . Initially, each cluster is a singleton composed of an individual cell in the cover. Then, in Line 4, we select the cluster which has the smallest (noisy) total count, and merge with its closest neighboring cluster in Line 5-8. The distance between two clusters are measured with the physical distance between their cluster centers 5 . In each step, the noisy total count of each cluster is computed as the sum of all (noisy) total counts of each cell within the cluster (Line 7). Since the total counts of cells are noisy, Algorithm 3 preserves DP.

Perturbing aggregated time series

To perturb aggregated time series, we build on the Fourier Perturbation Algorithm (FPA) [START_REF] Rastogi | Differentially private aggregation of distributed time-series with transformation and encryption[END_REF]: 5 The cluster center is the centroid of its consitituent cell polygons.

Algorithm 3 Cluster cells

Input: Minimum cover C = {c 0 , c 1 , . . . , c |C| }, Minimum total count τ , Noisy total counts Ŝ of cells in C Output: Cluster configuration E ⊂ 2 |C| 1: E := {{c 0 }, {c 1 }, {c 2 }, . . . ,

Pad

Fk to be a n-dimensional vector by appending n -k zeros, which is denoted by PAD n ( Fk ). Finally, the inverse DFT is applied to PAD n ( Fk ) to obtain a noisy version of X.

FPA provably guarantees ε-DP [START_REF] Rastogi | Differentially private aggregation of distributed time-series with transformation and encryption[END_REF]. Enhanced FPA [START_REF] Ács | Differentially private histogram publishing through lossy compression[END_REF] improves basic FPA by selecting the coefficients to be removed more effectively. Specifically, in Step 2, EFPA chooses k probabilistically using the exponential mechanism such that the values of k which minimize the root-sumsquared error E||X -X||2 ≤ n i=k+1 |Fi-1| 2 + 2k∆ 2 (X) ε (RSSE) have exponentially larger probability to be selected. In this paper, we improve the accuracy of EFPA in two ways. First, instead of the Discrete Fourier Transform (DFT), we Algorithm 4 EFPAG Input: Truncated time series X with length n, Privacy budget ε, δ, L 2 -sensitivity of X: ∆ 2 (X) Output: Noisy time series X with length n

1: F := DCT(X) // Discrete Cosine Transform 2: Compute u G (X, k) = n i=k+1 |F i-1 | 2 + √ 2∆ 2 (X) √ k ln 1/2 (4/δ) ε for all 1 ≤ k ≤ n 3: Select k with probability ∝ exp -ε•u G (X,k) 4∆ 2 (X)
4: apply the orthonormal version of the Discrete Cosine Transform (DCT), which tends to provide smaller high frequency components due to its different boundary conditions [START_REF] Strang | The discrete cosine transform[END_REF]. This can result in smaller RSSE when these components are removed in Step 2. Moreover, since we use orthonormal DCT, the resultant scheme also preserves ε-DP [START_REF] Ács | Differentially private histogram publishing through lossy compression[END_REF]. Second, instead of adding Laplace noise, we add properly calibrated Gaussian noise to the first k Fourier coefficients of X thereby providing larger accuracy at the cost of weaker privacy. More specifically, when Fk is generated in Step 3, we employ the Gaussian mechanism instead of LPA and add i.i.d Gaussian noise G(σ) to each coefficient in F k , where σ = √ 2∆2(X) ln 1/2 (4/δ)/ε to provide (ε, δ)probabilistic DP based on Theorem 3. In addition, we select k with probability ∝ exp

Fk := F k + G( √ 2∆ 2 (X) ln 1/2 (4/δ)/ε) k 5: return X = IDCT(PAD n ( Fk )) //
-ε•u G (X,k) 4 in Step 2, where uG(X, k) = n i=k+1 |Fi-1| 2 + √ 2∆ 2 (X) √ k ln 1/2 (4/δ) ε
which follows from Theorem 5 in [START_REF] Ács | Differentially private histogram publishing through lossy compression[END_REF] and Theorem 3. When we use GPA in Step 3, the new scheme (with DCT) is denoted by EFPAG in the rest.

Since Gaussian noise has smaller variance than Laplacian, EFPAG provides better accuracy than EFPA. Specifically, the variance of the Gaussian noise added to the Fourier coefficients is 8∆2(X) 2 ln(4/δ)/ε 2 , which is independent of the number of retained coefficients (i.e., k). By contrast, the variance of the Laplace noise added to the coefficients in EFPA is 8∆2(X) 2 k/ε 2 which linearly increases with the number of retained coefficients. However, this improvement also leads to some privacy degradation which is measured by δ. In the sequel, we fix δ to 2 • 10 -6 (≤ 1/N ).

EFPAG is summarized in Algorithm 4, where the total budget ε is uniformly divided between GPA (Line 4) and exponential mechanism (Line 2), therefore, EFPAG is (ε, δ)probabilistically DP due to Theorem 1. Finally, in order to employ EFPA(G), we need to compute the L2-sensitivity of the counts in the cover C, i.e., ∆2(X C ).

Indeed, since E is a partitioning of C, ∆2(X E ) = ∆2(X C ).
Recall that, as a result of pre-sampling, at most a single visit of a user is retained in any slot, and at most ℓ visits per user over the whole week. This means that, for any t, there is only a single tower whose count can change (by at most 1) by modifying a single user's data. From Formula (1), it follows that the total change of all IRIS cell counts is at most 1 at any t, and hence ∆2(X 3 illustrates the improvement of our approach (Clustering + EFPAG) over simple LPA.

L ) ≤ ∆2(X T ) = √ ℓ based on the definition of L2-norm. Since C ⊆ T ∪ L, ∆2(X C ) ≤ √ ℓ. Figure

Further improvements

Although our approach is clearly superior to LPA, Figure 4 still suggests a large error on average. This difference between X and X is the result of two errors: the sampling error (between X and X) is attributed to pre-sampling, whereas the perturbation error (between X and X) is due to our perturbation scheme.

As illustrated by Figure 4a, sampling error mainly distorts large counts: although the noisy counts are close to the counts of the truncated (pre-sampled) time series between 9:00 AM and 11:00 PM, it is still far from the original count values. This significantly increases MRE.

In addition, as Figure 4b also shows, noisy counts also deviate from pre-sampled as well as from original counts around the local minimas (close to 4:00 AM every day), which further deteriorates MRE. This perturbation error is caused by the higher frequency components that are retained and perturbed by EFPA(G).

To alleviate these errors, we propose two further improvements: first, we improve the perturbation of total cell counts (Line 1 in Algorithm 2), which is used in cell clustering (Algorithm 3) and scaling (Line 6 in Algorithm 2). Then, as a post-processing step, we smooth out small counts (i.e., between 0:00 and 6:00 AM) through non-linear least-square fitting to diminish perturbation error.

Improved scaling

Recall that we scale back the normalized aggregated time series with Ŝi in Line 6 (Algorithm 2), where Ŝi = 167 t=0 X i t + L(2ℓ/ε). Since X i is the pre-sampled time series of cell i, Xi (Line 6) will be a scaled down version of the original time series X i due to the fact that the ℓ visits per individual are sampled uniformly at random. Also, as we have discussed in Section 5.1, adding Laplace noise directly to the original total count 167 t=0 X i t is very inaccurate, as the sensitivity of 167 t=0 X i t is ∆1(X) which is too large. We rather perturb the original total count 167 t=0 X i t using a different approach: we first approximate the relative frequencies of each tower by another constraint sampling, and scale back these frequencies to count values with the (noisy) total number of visits in the dataset. The main idea is that sampling requires only a small amount of noise to guarantee privacy, while the total number of all visits is so large that it tolerates a large noise magnitude.

In particular, we estimate the histogram H where a bin Hj represents the frequency of visits at tower j, i.e., Hj = 167 t=0 X j t /K, where K denotes the total number of tower visits in the dataset (K = 167 t=0 T ∈T X T t ). To do so, we sample a single visit per user uniformly at random, and create a new histogram H from the sampled visits (with size N ). Using this approximative histogram H, the total number of visits 167 t=0 X j t of a tower j is computed as ( Hj + L(2/ε ′ )) × K, where K = K + L(2∆1(X)/ε ′ ) and ∆1(X) is universally fixed for all input dataset 6 . Finally, having the noisy 167 t=0 X j t for each tower cell j, we can also compute the noisy 167 t=0 X L t for any IRIS cell L (using Formula (1)) and calculate Ŝi for all cell i in C. This technique is 2 • (ε ′ /2) = ε ′ differential private based on Theorem 1.

Figure 6a compares the accuracy of our sampling approach to perturb the relative frequencies of each tower (i.e., sampling is followed by adding L(2/ε ′ ) to each Hj) with the naive Laplace approach (i.e., L(2∆1(X)/ε ′ ) is added to each Hj without sampling). Sampling clearly boosts the accuracy of histogram perturbation (Figure 6a) especially for smaller values of ε ′ , and eventually yields significantly more accurate estimation of 167 t=0 X T t for all towers T (Figure 6b). The effect of scaling is illustrated in Figure 5a. Recall that the full privacy budget ε is divided equally between EFPA(G) and scaling (see Algorithm 1). Hence, our sanitization scheme is ε-DP (or (ε, δ)-prob. DP with EFPAG) based on Theorem 1.

Smoothing

In order to smooth out low (noisy) counts around the local minimas (around 4:00 AM each day), we fit an exponential curve to the noisy counts between 0:00 AM and 4:00 AM 6 ∆1(X) is fixed to 732 in this paper. Notice that although ∆1(X) is large so is K: K = 137, 255, 052 in our dataset, and |K -K|/K is less than 10 -5 on average where the counts are exponentially decreasing, and another exponential curve between 4:00 AM and 6:00 AM, where the counts are exponentially increasing. In particular, we fit function g(x, a, b) = a•exp(b•x) to the noisy counts, i.e., compute parameters a, b such that the error i ( Xi -g(xi, a, b)) 2 is minimized where xi runs over the hours of the given time intervals for each day, and then replace the noisy counts with the values of the fitted function. This is a standard non-linear least square fitting problem which can be approximated with any numerical minimization method (e.g., Levenberg-Marquardt algorithm [START_REF] Levenberg | A method for the solution of certain non-linear problems in least squares[END_REF]). Since this operation is performed on the noisy time series, it is already private. The effect of smoothing in illustrated in Figure 5b.

PERFORMANCE EVALUATION

We evaluate the utility of our scheme depending on the guaranteed privacy (i.e., ε) with EFPA and EFPAG, where δ = 2 • 10 -6 < 1/N . The minimum total count τ used in Algorithm 3 is adjusted such that the expected RSSE is less than 1% of the total count when all coefficients are retained in EFPA(G). That is, τ = √ 168 • σ/0.01, where σ 2 is the variance of noise added to each coefficient. We compare our approaches to the naive Laplace mechanism (LPA) that adds L(ℓ/ε) noise to each count of each time series in cover C. We use the CDR dataset described in Section 4.

First, we analyze the utility depending on the presampling size ℓ. Then, we show how pre-sampling combined with the improved scaling and smoothing boost accuracy, and also report the error distribution among individual IRIS cells. Finally, we measure the Earth Mover's Distance (EMD) which captures the error between spatial distributions in terms of geographical distances.

Utility depending on the pre-sample size

Recall that the number of visits retained per user (i.e., ℓ) determines the injected noise in the perturbation phase. In general, larger values of ℓ imply larger noise, which degrades utility. On the other hand, smaller values of ℓ preserve more information about individuals which results in a more accurate representation of the original dataset. The goal is to select a value of ℓ which yields the best trade-off. Nevertheless, we experimentally show next that our scheme exhibits stable performance with quite different values of ℓ.

In Figure 7, we report the utility of our scheme with EF-PAG for three values of ℓ: 10, 30 and 168. In each case, perturbation is followed by the improvements described in Section 5.4. Specifically, we computed the average of MRE and Pearson correlation over all cells. We repeated the whole process 20 times and plotted the mean and standard deviation of the average MRE and PC over all executions. LPA has significantly larger error for all values of ℓ, and provides especially poor results for smaller values of ε. By contrast, our scheme does not only provide practical utility even for stringent privacy guarantee, but also has stable performance for different ℓ. For instance, for ε = 0.3, MRE is less than 20%, while PC is larger than 0.95. Therefore, the output of our scheme has almost perfect linear correlation with the original data thanks to the combination of clustering and the Fourier perturbation approach. Moreover, the variations of these values are very moderate: 0.2 ± 0.03 and 0.95 ± 0.03, respectively, for different values of ℓ. In the rest of the paper, we fix ℓ to 30.

EFPAG and EFPA 7 are compared in Figure 8. EFPAG outperforms EFPA especially for smaller ε: MRE is reduced by 0.03 and PC is increased by 0.02 on average.

Pre-sampling with scaling

The distortion effect of pre-sampling is mitigated by the improved scaling step detailed in Section 5.4.1. Our aim now is to show that scaling and smoothing indeed results in better utility. Figure 8 depicts a variation of our scheme, denoted by "EFPA-" when the improvements described in Section 5.4 are not employed after perturbation. The results show that MRE is reduced by 0.07 on average when the presample size is diminished to ℓ = 30 and improvements are employed. By contrast, PC is increased only by about 0.01 for smaller values of ε; the change is not so significant due to the fact that scaling does not influence linear correlation, 7 Recall that EFPAG adds Gaussian noise whereas EFPA adds Laplacian noise to the retained Fourier coefficients. and smoothing modifies relatively small number of counts in general.

In Figure 9, we also plotted the average relative error depending on the time for our scheme with and without improvements. In particular, we computed the relative error and took the average over all cells in each hour. Figure 9a confirms that scaling significantly diminishes the relative error in daylight when counts are larger. The improvement can be almost a factor of 4. This has particular importance in practice, as location counts in daylight are usually more important than at night. The MRE and PC of individual IRIS cells are also illustrated by color maps in Figure 11. This figure shows that our scheme can provide practical utility for most cells with strong privacy guarantee. Specifically, the average MRE over all cells is only 0.17 with ε = 0.3.

Error variation among IRIS cells

Earth Mover's Distance

In order to compare the sanitized spatial probability distribution with the original one at a given time, we use Earth Mover's Distance (EMD) [START_REF] Rubner | The earth mover's distance as a metric for image retrieval[END_REF]. EMD measures the "amount of energy" (or cost) needed to transform one distribution to another, and is a metric for probability distributions (i.e., location counts have to be normalized). Formally, for any t, EMD(X L t , XL t ) = min {f ij } i,j fijdij such that fij ≥ 0,

j fij ≤ X i t / k X k t , i fij ≤ Xj t / k Xk t
, where {fij} denotes the set of all possible flows (each fij represents the amount of probability mass transported from IRIS cell i to j), and dij is the geographical distance between the centers of cells i and j, resp. Intuitively, EMD measures the meters of error between two spatial density maps. Figure 9b reports the EMD depending on the time. The mean EMD over the whole week is 258 meters for EFPA, 188 meters for EFPAG, and 341 meters for LPA.

CONCLUSIONS

The goal of this work is to demonstrate through a realworld application that differential privacy can be a practical model for data anonymization, even if the input dataset has large dimension and/or is highly sensitive. We showed that, in order to achieve meaningful accuracy, the sanitization process has to be carefully customized to the application and public characteristics of the dataset. We strongly believe that there are no "universal" sanitization solutions that fit all applications, i.e., provide good accuracy in all scenarios. In particular, achieving the best performance requires to find the most faithful and concise representation of the data, such that it withstands perturbation. In our application (i.e., spatio-temporal density), clustering and sampling with Fourier-based perturbation are seemingly the best choices due to the periodic nature and large sensitivity of location counts. We experimentally showed that our scheme can provide practical utility and strong privacy guarantee.

Theorem 4 (

 4 [START_REF] Mcsherry | Mechanism design via differential privacy[END_REF]) Given a utility function u : (D × R) → R for a database D, the mechanism A, A(D, u) = return r with probability ∝ exp εu(D, r) 2∆ugives ε-differential privacy, where ∆u = max ∀r,D 1 ,D 2 |u(D1, r) -u(D2, r)|.

Figure 1 : 3 :

 13 Figure 1: IRIS cells of Paris (left) and Voronoi-tesselation of tower cells (right)

  in O(|E| |V |), where |V | = |T| + |L| and |E| ≤ |T||L|.

Figure 2 : 2

 22 Figure2: Covering Paris with the largest tower and IRIS cells. Each cell is colored based on their mean count. Large cells tend to have large counts (closer to red), while small cells are less populated (closer to blue). The minimum cover (Figure2c) includes large, more populated cells: IRIS cells from the city center, and tower cells around the perimeter.

Figure 3 :

 3 Figure 3: Noisy time series of an IRIS cell (ε = 0.3, ℓ = 30)

Figure 4 :Figure 5 :

 45 Figure 4: Our scheme before improvements (ε = 0.3, δ = 2 • 10 -6 , ℓ = 30).

  Perturbing H with sampling and without sampling.

  MRE of the noisy 167 t=0 X T t over all towers T .

Figure 6 :

 6 Figure 6: Perturbing the total visits 167 t=0 X T t of each tower cell T .

Figure 7 :Figure 8 :

 78 Figure 7: Utility depending on the pre-sampling size ℓ.

Figure 9 :Figure 10 :

 910 Figure 9: Error depending on time with EFPAG (ε = 0.3, ℓ = 30)

Figure 10

 10 Figure 10 shows through box plots how MRE and Pearson correlation change among IRIS cells; we compute these metrics for each IRIS cell, and compute the corresponding box plot over the metric values of all cells. Although medians do not change significantly for different values of ε, MRE has larger variation for smaller ε, i.e., there are more cells which have larger error.The MRE and PC of individual IRIS cells are also illustrated by color maps in Figure11. This figure shows that our scheme can provide practical utility for most cells with strong privacy guarantee. Specifically, the average MRE over all cells is only 0.17 with ε = 0.3.

Figure 11 :

 11 Figure 11: MRE and PC of each IRIS cell (ℓ = 30, ε = 0.3)

  {c |C| }} 2: ŜE := c i ∈E Ŝi for each cluster E ∈ E 3: while ∃E ∈ E such that ŜE < τ do 4: E := arg min E∈E ŜE 5: Let E ′ be the geographically closest neighbor of E 6: E := E ∪ E ′ 7:

	ŜE := ŜE + ŜE ′ 8: Remove E ′ from E 9: end while
	1. Compute	the	Fourier	coefficients	F	=
	F0, F1, . . . , Fn-1 of the input aggregated time
	series X with length n by discrete Fourier transform.
	3. Generate the noisy version of F k , denoted by Laplace mechanism: add i.i.d Laplace noise L( Fk , by √ k/ε) to each coefficient in F k .

2. Remove the last n -k coefficients from F, which correspond to the high-frequency components in X, and retain only the first k elements of F, denoted by F k . Note that k is an input to the algorithm.

  Inverse DCT

		350
		300	Original Private (MRE: 0.73, PC: 0.59)	250	Original Private (MRE: 0.38, PC: 0.98)
		250		200
	Visit count	150 200	Visit count	100 150
		100
		50		50
		M o 4 12 20 4 12 20 4 12 20 4 12 20 4 12 20 4 12 20 4 12 20 T u W ed T hu Fr i Sa t Su n 0	M o 4 12 20 4 12 20 4 12 20 4 12 20 4 12 20 4 12 20 4 12 20 T u W ed T hu Fr i Sa t Su n 0
			(a) LPA		(b) Our approach

National Institute of Statistics and Economics: http://www.insee.fr/fr/methodes/default.asp?page= zonages/iris.htm

Available on IGN's website (National Geographic Institute): http://professionnels.ign.fr/contoursiris

http://www.orange.com

Tower cells are represented by their Voronoi polygons as it is depicted in Figure1