Erwan Bousse
email: erwan.bousse@gmail.com

David Mentré
email: d.mentre@fr.merce.mee.com

Benoît Combemale
email: bcombema@irisa.fr

Benoît Baudry
email: benoit.baudry@inria.fr

Takaya Katsuragi
email: katsuragi.takaya@ap.mitsubishielectric.co.jp

Aligning SysML with the B Method to Provide V&V for Systems Engineering

Systems engineering, and especially the modeling of safety critical systems, needs proper means for early Validation and Verification (V&V) to detect critical issues as soon as possible. The objective of our work is to identify a verifiable subset of SysML that is usable by system engineers, while still amenable to automatic transformation towards formal verification tools. As we are interested in proving safety properties expressed using invariants on states, we consider the B method for this purpose. Our approach consists in an alignment of SysML concepts with an identified subset of the B method, using semantic similarities between both languages. We define a restricted SysML extended by a lightweight profile and a transformation towards the B method for V&V purposes. The obtained process is applied to a simplified concrete case study from the railway industry: a SysML model is designed with safety properties, then automatically transformed into B, and finally imported into Atelier-B for automated proof of the properties.

INTRODUCTION

Systems engineering is a field that involves many different stakeholders, each of them with a different role in the development and with a specific set of skills. This implies handling different views of the system, each of them adapted to the specialist it is dedicated to. For that matter, Domain Specific Modeling Languages (DSML) offer a satisfying solution: designing languages whose syntax and purposes exactly fit stakeholder roles, while being adapted to their domains. For instance, the Systems Modeling Language (SysML) is a graphical language adapted to systems engineering that can be used throughout a system development. A particular kind of systems are safety critical systems, for which failures can lead to casualties or the loss of equipment.

Validation and Verification (V&V) of safety critical systems is required, as safety properties and compliance with standards such as EN 50128 (railway applications), IEC 61508 (electronic safety related systems), or DO 178C (softwarebased aerospace systems) must be ensured. This must be done as early as possible during the design phase to limit rollbacks and underlying costs. In this regard, formal methods can be used to check specified properties of the system in all possible configurations. However, existing approaches require the use of formal notations that can be difficult to read and understand, especially for experts of a particular field. A way to provide V&V for existing DSMLs is to provide a translation towards formal notations, so that existing theorem provers or model checkers can be used.

Our work focuses on what we call the alignment of a DSML with a formal language. The idea is to identify good constructs and good practices of the formal language, then to search semantic similarities between such constructs and the considered DSML, and finally to define a transformation from the DSML to the formal language using Model Driven Engineering (MDE) techniques. We apply these principles to SysML -using the Object Constraint Language (OCL) for safety properties and the Action Language for Foundational UML1 (Alf) for behaviors -and the B method. The resulting process relies on SysML block definition and state machine diagrams on the one hand, and on B developed modules linked by imports and sees links on the other hand. This process leads to concrete proved software using the Atelier B tool. To complete this work, we implemented the first version of a prototype in Kermeta that realizes this transformation, and we provide a concrete application of this tool-supported process to a simplified example of railway specification provided by Mitsubishi Electric.

Remaining sections are organized as follows. Section 2 lays down the industrial context of this work at Mitsubishi Electric. Section 3 introduces languages we consider in our approach -namely SysML and the B method. Section 4 presents our contribution, which consists in an alignment of SysML with the B method. Sections 5 and 6 conclude with related work and final words.

INDUSTRIAL CONTEXT

In this second section, we present the context of this work at Mitsubishi Electric and the underlying objectives.

Goals and technical choices

The scope of this work at Mitsubishi Electric is the design of safety-critical railway systems; in particular, we have the following goals. (1) Modeling languages must be adapted to systems engineering and its practitioners, which may not be formal methods specialists. (2) We want to be able to design systems at least at the implementation level, i.e. models translatable into concrete deployable embedded software.

(3) It is necessary to manage good traceability between informal requirements and the modeled system, especially for safety properties. (4) Ultimately, we must be able to prove such properties using formal methods with tool support, keeping in mind scalability issues. ( 5) Finally, we need to consider at least safety properties expressed as invariants on states of the system. The SysML language [START_REF] Omg | OMG Systems Modeling Language (OMG SysML)[END_REF] is a rising modeling language for systems with dedicated traceability specific features: a requirements diagram and links between requirements and elements of the model. Additionally, this language is complete enough to provide ways to design systems at any level of abstraction, including the implementation level. Following goals 1, 2 and 3, we choose to use SysML for the modeling part of our process.

As the B method [START_REF] Katsuragi | The B-book: assigning programs to meanings[END_REF] has been used in important safetycritical system realizations in Europe, such as the Paris Métro Line 14 [START_REF] Behm | Météor: A successful application of B in a large project[END_REF], we looked at this formal method with interest. The B method relies on theorem proving, which means that scalability is rather well handled, unlike state explosion problems that may appear with model checking techniques [START_REF] Demri | A parametric analysis of the state-explosion problem in model checking[END_REF][START_REF] Clarke | My 27-year quest to overcome the state explosion problem[END_REF]. In terms of granularity, this method goes down to the implementation of the system. B notations and concepts can be difficult to apprehend, but this is not a problem since we only want to use this formal method for V&V purposes. Overall, it meets satisfactorily goals 2, 4 and 5

Case study: a railway crossing controller

For our work, we use a proprietary case study provided by Mitsubishi Electric. It is a technical specification and requirements document describing a railway crossing controller. Figure 1 shows the configuration of the crossing: there are two lines (inbound and outbound), both with a critical section between two train sensors, and cars may pass on the crossing road. Each track has two sensors: one called ADC that detects trains arriving in the warning section, another called BDC for trains leaving the warning section. The main requirement is to detect when trains are in the critical sections in order to activate the barriers, the bells, and the lamps accordingly. Being a life-critical system, properties such as "when there is a train in a critical section, lamps are lit" or "barriers are closed or closing when there are trains in the critical sections" have to be proved.

BACKGROUND 3.1 SysML, OCL and Alf

The Systems Modeling Language (SysML) [START_REF] Omg | OMG Systems Modeling Language (OMG SysML)[END_REF] is a graphical modeling language managed by the Object Management Group (OMG). It is oriented towards general purpose systems engineering applications. It is well suited for systems mixing hardware and software with constraints of many kinds, as in rail systems. Over the 9 SysML diagrams, we only use in this paper the block definition diagram (structure modeling) and the state machine diagram (reactive behavior modeling), a block being the basic building block of the language to define system elements. It is also possible to design extensions for SysML using its profile mechanism.

Additionally, we use the possibility provided by the SysML language to enrich models using other languages. The Object Constraint Language (OCL) is a specification language that allows to enrich models with side-effect-free annotation; we use it to write invariant constraints in our SysML models. The Action Language for Foundational UML (Alf) [START_REF] Omg | Action Language for Foundational UML (Alf), Concrete Syntax for a UML Action Language, FTF -Beta 1[END_REF] is an action language with a Java-like syntax for describing behaviors, such as values assignments or operation calls; we use it to define low level behaviors in our SysML models.

The B Method

The B method [START_REF] Katsuragi | The B-book: assigning programs to meanings[END_REF] is a field-proven formal method that was successful in various safety-critical rail system applications like Paris Métro Line 14 or Roissy VAL [START_REF] Behm | Météor: A successful application of B in a large project[END_REF][START_REF] Badeau | Using B as a high level programming language in an industrial project: roissy VAL[END_REF][START_REF]Formal methods in industry: achievements, problems, future[END_REF]. Tools are available for operational use like Atelier B [START_REF] Clearsy | B Language Reference Manual -Version 1[END_REF], a proprietary integrated development environment developed by ClearSy. For our work, we consider Atelier B as a part of the suggested process: we generate B code with Atelier B specific syntax.

A B project is a collection of modules that models a system. A module is a part of the system, and consists in a stack of components which can be abstract machines, refinements or implementations. An abstract machine is the visible specification of a module. A refinement refines either an abstract machine or another refinement, and there can be several refinement steps. An implementation is the final refinement step of a module. Most of the B language is a specification language, which includes set theory, first order logic, indeterminism and abstract data. This specification language is used to define abstract machines and refinements. Another part of the language, called B0, is a formally defined programming language used in implementations that can be compiled to C or Ada. Each module is made of at least one abstract machine, zero or more refinements, and zero or one implementation. A module only made of an abstract machine is said to be abstract, while a module with an implementation is developed. It is possible to instantiate modules, which allows the use of instances of lower level modules within the implementation of a module. This is done by using the imports link. On the other hand, the includes link is a way to decompose an abstract machine into a set of abstract machines, but without instantiations: this is equivalent to splitting the machine into multiple files instead of a single one. Finally, the sees link is a read-only access from a module to another; it is required when data has to be shared between modules, e.g. constants definition.

In this paper, even though B has abstraction capabilities, we consider B models that describe implementations of systems. Specifications of such systems are translated in two kind of B elements: invariant properties to be proved on the implementations, and B abstractions corresponding to assumptions on the environment.

ALIGNING SYSML WITH B

In this section we present our main contribution: based on a B subset we first delimit, we identify a subset of SysML aligned with this subset of B, and we define of a mapping between the two languages we eventually obtain.

Industrial use of the B method

There always are multiple ways to model a system with a specific language; however, all languages have what we can call "good practices", which are rules to apply in order to obtain "good" models or programs. For instance, in objectoriented programming, design patterns give solutions considered effective for some precise problems, like strategy for algorithm modularity. In the case of formal languages, another factor than effectiveness comes into account: the ability to prove required properties of the model. The more constrained, simple and small a model is, the easier it is to analyze it. Applying these criteria can limit the state-space size of the model or give clearer relationships between model elements. Formal models must also be well decomposed, so that properties that are specific to a part of a system are easier to prove and do not depend on unrelated parts.

Regarding the B method, there mainly are two aspects to consider: what kinds of B components and what kinds of links between these components to use. We deduct from past industrial uses [START_REF] Badeau | Using B as a high level programming language in an industrial project: roissy VAL[END_REF] of the B method that the most important is to have well decomposed B projects: for scalability concerns, a property to prove must be as much as possible contained into a subpart of the model. The imports link is a way to do such a modular design, by defining a module using lower level instances of other modules. This implies that we work with developed modules instead of purely abstract modules. Refinement components are rarely used 2 : each abstract machine is refined only once into an implementation. Additionally, sees links are required when data has to be shared between modules, like global constants. To sum up, we restrict our usage of the B method to developed modules with one abstract machine and one implementation linked by imports and sees links. Concerning data, we only use primitive types (boolean, integer and enumerations) in this first approach.

Semantic similarities

2 Refinements are only used when an algorithm or a data structure are too complex to be proved directly at implementation level. After studying and restricting the B method, the idea is now to align SysML with this delimited subset. We call semantic similarities between languages features of both languages that are close semantically. In a nutshell, we want to be able to construct models in SysML that can be translated in this subset of the B method while preserving semantics.

Parallel read of the languages specifications

Searching for semantic similarities consists in reading official specifications of both languages and highlighting potentially related parts. It is important to note that the role of each concept is defined using natural language, which means that the parallels depend on our reading and understanding. This must not be confused with the fact that concepts may be defined using formalisms, like the B language with the B-Book [START_REF] Katsuragi | The B-book: assigning programs to meanings[END_REF], while we look at the purpose of the concepts.

To illustrate our approach, here is how we proceed for the first three similarities. Table 1 shows quotes of SysML and B specifications, each line putting side by side similar concepts of SysML and the B method. We find the following parallels.

A B project has the same purpose as a set of SysML blocks, which is the definition of a system. Such a block is very similar to a B module to model an element of the system. Finally, an imports link defines a subsystem of a module the same way a SysML part property does for a block.

Altogether, we identify 12 semantic similarities between SysML and the B method for this first approach -see Table 2 for the complete list. The complete analysis is available in [START_REF] Bousse | Requirements management led by formal verification[END_REF]. Using these similarities, we define a subset of SysML which allows only a few constructs among the many possibilities of SysML. In short, each model is a collection of blocks with primitive data in value properties, linked by part property links for the system decomposition, with textual operations to define their behavior and required constraints over variables. Such models can be entirely designed with block definition diagrams.

Additionally, we use opaque behavior and opaque expression features of SysML to enrich models with textual behaviors and expressions written in Alf and OCL respectively. This part of the translation relies on straightforward parallels between B operators or instructions and Alf, and between B predicates and OCL. We do not present these semantic similarities in this paper (refer to [START_REF] Bousse | Requirements management led by formal verification[END_REF] for more details). "Blocks provide a general-purpose capability to model systems as trees of modular components. The specific kinds of components, the kinds of connections between them, and the way these elements combine to define the total system can all be selected according to the goals of a particular system model" "A B module models a sub-system; it forms a part of a B project." "A Block is a modular unit that describes the structure of a system or element. " "Import is used to structure a B project into layers, since the implementation of a module is implemented by importing other modules."

"SysML blocks [. . . ] provide the ability to represent a system hierarchy, in which a system at one level is composed of systems at a more basic level. [. . . ] A part property holds instances that belong to a larger whole."

Some additions to our subset of SysML

Even with a subset of the B method and a subset of SysML that are semantically similar, we are not completely satisfied yet for two reasons. First, we did not find SysML counterparts for some essential concepts of the B method: the notion of main module of a project, being able to define abstract elements, and how to consider invariants defined in the implementation of a module. Second, we want to be able to design reactive systems in a practical way.

To overcome the first reason, we use UML profile mechanism to define an very simple extension to the SysML language. Firstly, a main stereotype to tag the main block of the SysML model. Secondly, an abstract stereotype to tag value properties that are abstract and should not appear in the final implementation -this is necessary to model hypotheses on the environment. Thirdly, a gluing stereotype to tag constraints that both affect properties of a block and properties of its parts -this corresponds to B invariants expressed in the implementations of modules.

For the second reason, we rely on the SysML state machine feature to define reactive behaviors of blocks. Since there is no such concept within the B method, this part does not rely on any semantic similarities. However, preserving semantics in a translation is nonetheless possible, as shown by Sekerenski [START_REF] Sekerinski | Graphical design of reactive systems[END_REF] who defined rules to transform UML state machines into B abstract machines. Using and adapting this work, we add a subpart of SysML state machine concepts to our restricted SysML -putting aside most pseudostates and entry/do/exit behaviors amongst other things.

We defined our customized SysML with14 rules (presented in [START_REF] Bousse | Requirements management led by formal verification[END_REF]) written in natural language that restrict SysML usage. Further work may lead us to write them using the OCL language over the SysML metamodel.

Modeling the case study

To experiment with the restricted SysML we eventually define, we model a simplified version of the case study introduced Section 2.2. We consider that only two lamps have to be lit when trains are in the warning section between the sensors and that a failure lamp must be activated when a failure state is reached. troller, whose properties are a counter train_counter for trains in the warning section and a constant MAX_TRAINS for the maximum number of trains that can exist in the warning section at a given time. Two operations are defined: ADC() called when a train arrives in a warning section and BDC() when a train leaves one. The Lamp block models a physical lamp: it has a boolean is_on set to true (respectively false) when a turn_on() (respectively turn_off()) operation is called. This variable is abstract because we do not need at the implementation level to keep in memory the state of the lamp: it is an hypothesis we consider on the behavior of a real lamp. Operations turn_on/turn_off concretely act on these lamps using a hardware API, which is not shown here. The Controller block possesses three part properties lamp1, lamp2 and failure_lamp, which are instances of the Lamp block. Initial conditions are the following: no trains on the tracks, lamps are not lit. Two safety properties are defined within the Controller block using gluing constraints: the first one guarantees that lamps are lit when the train counter is strictly positive; the second one ensures that the failure lamp is lit when the train counter becomes inconsistent.

A state machine diagram representing the reactive behavior of the Controller block is also part of this SysML model, but due to the page limit, we chose in this paper to focus on the part of our SysML subset based on semantic similarities. Refer to [START_REF] Bousse | Requirements management led by formal verification[END_REF] for the complete case study.

Translation of SysML into B

Now that we have a SysML subset corresponding to the targeted B subset, the last step of our work is the definition of a translation toward the B method. Semantic similarities we identified between both languages allowed us to delimit the SysML we use, but also directly gives us a mapping between SysML and B concepts, which directly leads to transformation rules. Figures 3 and4 are parts of a B model that results from the application of such rules of the SysML model introduced Section 4.2.3.

Each block of the SysML model is translated into a B module, which is made of an abstract machine with the same name as the block (ex. Lamp) and an implementation with a i suffix (ex. Lamp i). Part properties of the Controller block (ex. lamp1) are translated into imports links with the same names. Regular value properties (ex. train_counter) become B concrete variables, whereas those with their is-ReadOnly attribute set (ex. MAX_TRAINS) become B concrete constants. Abstract elements of the SysML model (is_on property, behaviors of the turn_on and turn_off operations) are translated into purely abstract B data and behaviors. Constraints of the SysML model are translated into their B predicate counterpart as a clause of a B invariant of the corresponding module. If the constraint is tagged with the gluing stereotype, we put it in the invariant of the implementation -this is necessary since B imports are only visible from there. Note that we do not present the part of the transformation that concerns state machines, as we chose to focus on SysML blocks in this paper. When importing the whole model in Atelier B, 49 Proof Obligations (PO) are generated 3 . They are all solved automatically. Our translation is defined by 19 rules established using minimal examples (presented in [START_REF] Bousse | Requirements management led by formal verification[END_REF]). Using Model Driven Engineering (MDE) techniques, we implemented a proprietary model transformation written in Kermeta [START_REF] Muller | Weaving executability into object-oriented meta-languages[END_REF]. The tool takes as input a SysML model encoded using the XML Metadata Interchange (XMI) format, and applies a modelto-model transformation defined between the SysML metamodel and the B metamodel 4 . Then, a model-to-text transformation is applied to the B model encoded in XMI to obtain a usable Atelier B project. More information about the tool can be found in a dedicated web page [START_REF] Bousse | SysML to B translator[END_REF].

RELATED WORK

Snook and Butler [START_REF] Snook | UML-B: Formal modelling and design aided by UML[END_REF] worked on a translation from UML towards the B method. While their motivation is similar to ours -extending UML to provide an easier way to use the B method-their solution is different. They transform a whole UML model into a single B module, while we prefer to have several machines within a well decomposed B project. To enrich their UML models, they designed a language called µB with a B-like syntax, while we chose to use existing OCL and Alf languages. Finally, they designed a UML profile to tag UML elements with many B concepts, while the profile we defined for SysML is as non-intrusive and lightweight as possible. We believe that our choices lead to a process requiring less knowledge of the B method, while taking into account how the B method is industrially used.

Many other works offered ways to transform UML models into B. Lano, Clark and Androutsopoulos [START_REF] Lano | Uml to b: Formal verification of object-oriented models[END_REF] worked on the translation of UML-RSDS models into B. Laleau and Polack [START_REF] Laleau | Coming and going from UML to B : a proposal to support traceability in rigorous IS development[END_REF] suggested a transformation process oriented towards Information Systems (IS) modeled in UML, using an extension of UML called IS UML. Meyer and Souquières [START_REF] Meyer | A Systematic Approach to Transform OMT Diagrams to a B Specification[END_REF] originally worked on the translation of Object Modeling Technique (OMT) models. However, all these approaches only worked with B abstract modules and includes links, while we focused and developed modules and imports links.

We already mentioned the work of Sekerinski on state machines [START_REF] Sekerinski | Graphical design of reactive systems[END_REF]. He is one of the few who went through most possibilities offered by state machines: composite states, orthogonality, spontaneous transitions, etc. However, he does not consider structure modeling with blocks or classes, and only provide a way to translate a state machine diagram into a single B abstract machine.

To our knowledge, the work of Laleau, Semmak, Matoussi and Gnaho [START_REF] Laleau | A first attempt to combine SysML requirements diagrams and B[END_REF][START_REF] Matoussi | A goal-based approach to guide the design of an abstract Event-B specification[END_REF][START_REF] Matoussi | An Event-B formalization of KAOS goal refinement patterns[END_REF] is the only other attempt to combine SysML and the B method in a same process. Their work focuses on the SysML requirement diagram and is based on the following observation: relationships between SysML requirements are quite limited compared to more elaborated goalbased approaches. Their solution consists in adding KAOS (Keep All Objectives Satisfied) concepts inside SysML by defining a SysML profile. Then, they define a transformation from this extended SysML towards Event-B -a recent event based alternative to the B method -which turns requirements into B events. While this work may lead to an interesting goal-based approach for SysML, it only focuses on SysML requirements and is quite far from our own objective of modeling complete functional systems using blocks and state machines.

CONCLUSION

In this paper, we looked at how to provide early V&V for the SysML language using the existing tool-supported B method. Our goals were to consider models designed at the implementation level and to use this formal method in a scalable way. Our contribution relies in the alignment of SysML with the B method in three steps. The first one is the definition of a subset of the B method that focuses on good decomposition, using imports and sees links between B modules. The second step is a search for semantic similarities between SysML and the previously identified subset of the B method, in order to define a subset of the SysML language semantically close to the subset of B. Our retained subset relies on block definition and state machine diagrams. We had to define a SysML profile in order to fill a gap with some B concepts. The third step relies on found similarities to define a transformation with a semantic gap as little as possible. We implemented it using the Kermeta language. We applied this SysML to B transformation on a concrete example from the rail industry, and managed to prove all safety properties.

Figure 1 :

 1 Figure 1: Case study: a railway crossing

Figure 2

 2 Figure 2 is a block definition diagram representing all blocks of the system and their links. The main block is the Con-

Figure 2 :

 2 Figure 2: Block definition diagram of the case study SysML model

Figure 3 :

 3 Figure 3: B module obtained from the Lamp block of the SysML model of Figure 2.

Figure 4 :

 4 Figure 4: B module obtained from the Controller block of the SysML model of Figure 2.

Table 2 :

 2 List of semantic similarities we consider

	B concept	SysML concept
	Project	Set of blocks
	Module	Block
	Imports link	Part property
	Sees link	Reference property
	Basic type	Value type
	Basic data	Value properties
	Constants	isReadOnly meta-attribute
	Enumerated set	Enumeration
	Component parameters	Default values
	Initialization and valuing Default values
	Operation	Operation (with behavior)
	Invariant	Constraint property

Table 1 :

 1 Semantic similarities between B method projects and SysML blocks to model structure Quote from ClearSy B Language Manual[START_REF] Clearsy | B Language Reference Manual -Version 1[END_REF] Quote from OMG SysML specification[START_REF] Omg | OMG Systems Modeling Language (OMG SysML)[END_REF] "A complete development in B corresponds to a B project.A project enables formally a system of any type. [. . . ] A B project refers to a complete set of B module instances. The components of these module instances are connected by links."

  REFINES Controller IMPORTS lamp1.Lamp, lamp2.Lamp, failure lamp.Lamp INVARIANT ((train counter > 0) ⇒(lamp1.is on=TRUE ∧ lamp2.is on=TRUE))∧ ((train counter < 0 ∨ train counter > MAX TRAINS) ⇒ failure lamp.is on=TRUE)

	MACHINE Controller
	SETS Controller states = {WSEmpty,WSnotEmpty,Failure} CONCRETE CONSTANTS	IMPLEMENTATION Controller i
	MAX TRAINS
	PROPERTIES
	MAX TRAINS ∈ INT ∧ MAX TRAINS = 15
	CONCRETE VARIABLES
	Controller state, train counter
	INVARIANT
	Controller state∈Controller states ∧ train counter∈INT∧ (Controller state = WSnotEmpty ⇒(train counter≤MAX TRAINS∧train counter>0))∧ (Controller state = WSEmpty ⇒ (train counter = 0)) INITIALISATION Controller state := WSEmpty train counter := 0 OPERATIONS ADC = ... BDC = ...	VALUES MAX TRAINS=15 INITIALISATION Controller state := WSEmpty ; train counter := 0 OPERATIONS ADC = ... BDC = ... END
	END

Unified Modeling Language

PO distribution: 2 for initial values, 12 for underflows/overflows of the counter, and 35 for safety invariants.

[START_REF] Behm | Météor: A successful application of B in a large project[END_REF] We wanted to use the B metamodel of Idani[START_REF] Idani | Dynamic graphical UML views from formal B specifications[END_REF], but as it focuses on dependencies between B constructs, we had to design one that considers the complete B abstract syntax.

This work is a first experiment limited to simple data structures and without making real use of the B method refinement capabilities. Aligning both languages aims to preserve semantics during the transformation, but such preservation is not proved thus far. Moreover, this work lacks effective support of traceability using SysML mechanisms and errors identified by the B method cannot be easily linked to the original SysML model. Further work will extend this V&V approach to other DSMLs and formal methods, with a particular focus on bidirectional transformation for traceability and reflecting identified errors on the original model.


Fabien Campillo

Frédéric Cérou

François Le Gland

Rivo Rakotozafy

Projet Mefisto

Particle and cell approximations for nonlinear filtering

Keywords: Nonlinear filtering, particle approximation, cell approximation, target motion analysis

We consider the nonlinear filtering problem for systems with noise-free state equation. First, we study a particle approximation of the a posteriori probability distribution, and we give an estimate of the approximation error. Then we show, and we illustrate with numerical examples, that this approximation can produce a non consistent estimation of the state of the system when the measurement noise tends to zero. Hence, we propose a histogram-like modification of the particle approximation, which is always consistent. Finally, we present an application to target motion analysis.

Approximations particulaire et cellulaire pour le filtrage non linéaire

Résumé : Nous considérons le problème de filtrage non-linéaire pour les systèmes sans bruit de dynamique. Nous étudions d'abord une approximation particulaire de la loi a posteriori, et nous donnons une estimation de l'erreur d'approximation. Nous mettons ensuite en évidence, et nous illustrons à l'aide d'exemples numériques, le fait que cette approximation peut donner un estimateur non-consistant de l'état du système, quand le bruit d'observation tend vers zéro. Nous proposons alors une modification de l'approximation particulaire, de type histogramme, qui est toujours consistante. Nous présentons enfin une application à la trajectographie passive.

Mots-clé :

Filtrage non-linéaire, approximation particulaire, approximation cellulaire, trajectographie passive.

Introduction

Consider the following nonlinear filtering problem Ẋt = b(X t ) , X 0 unknown,

z k = h(X t k ) + v k , (1)
where

t 0 < t 1 < • • • < t k < • •
• is a strictly increasing sequence of observation times, and {v k , k ≥ 0} is a Gaussian white noise with non singular covariance matrix R.

In this model, only the initial condition X 0 is unknown. The problem is to estimate X 0 , at any instant t k , given the past measurements z 1 , • • • , z k .

We can expect that, for this particular model, the nonlinear filtering problem will reduce to a parameter estimation problem of the unknown parameter X 0 . In this case, x 0 ∈ IR m will denote the true value of the parameter.

Whether we take a Bayesian approach or not, the goal is to compute :

Bayesian : the conditional probability distribution

µ k 0 (dx) = P (X 0 ∈ dx|Z k ) , (2)
where

Z k = σ(z 1 , • • • , z k ) , (3)
Non Bayesian : the likelihood function Ξ k (x) corresponding to the estimation of the unknown parameter X 0 ∈ IR m .

In this work, we focus on problems where the a priori available information on the initial condition X 0 is quite poor, for example X 0 ∈ K, where K is a compact subset of IR m . In this case, the two points of view -Bayesian and non Bayesian -are rather close.

We can study the asymptotic behavior (consistency, convergence rate, etc.) when (i ) the number k of observations tends to ∞ (long time asymptotics),

(ii ) or when the noise covariance matrix R tend to 0 (small noise asymptotics).

We also want to study the case where the system (1) is not identifiable. In this case the conditional probability distribution µ k 0 (dx) does not concentrate around the true value x 0 ∈ IR m , neither in the long time asymptotics, nor in the small noise asymptotics, but it concentrates around the subset M (x 0 ) ⊂ IR m of those points which are indistinguishable from the true value parameter. In this context, the relevant statistics that we should compute from the conditional probability distribution µ k 0 (dx) is not the usual pointwise Bayesian estimator of the conditional mean type :

X k 0 = x µ k 0 (dx) . We rather use the concept of confidence region : Definition 1.1 (Confidence region) For any given level 0 < α ≤ 1, a confidence region of level α is defined by :

D α k ∈ Arg min D∈D α k λ(D) , with D α k = {D ⊂ IR m : µ k 0 (D) ≥ α} , (4)
where λ denote the Lebesgue measure on IR m . This is an extension of the maximum a posteriori (MAP) estimator which can be seen as the limit, when the level α ↓ 0 , of the sequence of decreasing confidence regions. Our goal is to numerically compute the conditional probability distribution µ k 0 (dx).

In Section 2, we recall some properties of this nonlinear filtering problem, especially the crucial importance of the flow of diffeomorphisms {ξ s,t (•) , 0 ≤ s ≤ t} associated with the differential equation Ẋt = b(X t ), i.e.

X t = ξ s,t (X s ) , 0 ≤ s ≤ t , which allows an explicit formulation of the conditional probability distributions µ k 0 (dx) = P (X 0 ∈ dx | Z k ) and µ k (dx) = P (X t k ∈ dx | Z k ). As an example, we first consider the case where the initial probability distribution µ 0 (dx) = P (X 0 ∈ dx) is discrete, and the case where it is absolutely continuous with respect to the Lebesgue measure on IR m . Then we develop two numerical methods for the approximation of the conditional probability distribution µ k (dx). Section 3 is devoted to a particle-like approximation algorithm

µ k (dx) µ H k (dx) , with µ H k = i∈I a i k δ x i k ,
as a convex linear combination of Dirac measures, called particles. This kind of approximation has been introduced and studied by Raviart [START_REF] Raviart | An analysis of particle methods[END_REF] for the first order deterministic PDE's.

In the context of nonlinear filtering, this approximation technique has already been applied to various real case studies, see Campillo-Le Gland [START_REF] Campillo | Application du filtrage non-linéaire en trajectographie passive[END_REF], Le Gland-Pardoux [START_REF] Le Gland | Application du filtrage non linéaire oprtimal à la restitution d'orbite de transfert[END_REF] and Cérou-Rakotozafy [START_REF] Érou | Implémentation parallèle du filtrage non-linéaire pour la trajectographie passive[END_REF].

The problem is to define, at each time t k , the positions {x i k , i ∈ I} and the weights {a i k , i ∈ I} of the particles. The more natural choice is to define the approximation µ H k (dx) as the conditional probability distribution of X t k given Z k , with µ H 0 (dx) instead of µ 0 (dx) as an initial probability distribution. With this choice, the algorithm is the following

x i k = ξ t k-1 ,t k (x i k-1 ) and a i k = c k Ψ k (x i k ) a i k-1 , (5)
for all i ∈ I, where, by definition, Ψ k (x) is the likelihood function for the estimation, given the measurement z k , of the parameter

X t k ∈ IR m , that is Ψ k (x) = exp -1 2 z k -h(x) 2 R -1 with x 2 R -1 = x * R -1 x (6)
and c k is a normalizing factor. The only error is in the approximation µ 0 (dx) µ H 0 (dx), and it is then important to choose appropriately the initial points {x i 0 , i ∈ I} and weights {a i 0 , i ∈ I}. For this purpose, we use the approximation proposed in Florchinger-Le Gland [START_REF] Florchinger | Particle approximation for first order SPDE's[END_REF].

The main points of this section are the following :

(i ) From the numerical point of view, we obtain an estimate of the error

µ k (dx) -µ H k (dx) ,
in terms of the discretization parameter H, see Theorem 3.1 below.

(ii ) From the point-estimation point of view, the particle algorithm consists in restricting the parameter set to a set G H = {x i 0 , i ∈ I} ⊂ IR m of possible initial conditions : this is a misspecified estimation problem (i.e. the true value x 0 is not necessarily in G H ). In case where the model ( 1) is not identifiable the setwise Bayesian estimator (4) based on the particle approximation could be non consistent in both the long time asymptotics and the small noise asymptotics (for a fixed discretization parameter H). This phenomenon will be illustrated by simulation results.

Because of this possible non-consistency, we will study in Section 4 another approximation technique, where at time t k , instead of evaluating as in [START_REF] Érou | Implémentation parallèle du filtrage non-linéaire pour la trajectographie passive[END_REF] the value of the likelihood function at the point x i k , we evaluate a generalized likelihood function on a neighborhood B i k of the particle position x i k , see definition [START_REF] Kunita | Stochastic Flows and Stochastic Differential Equations[END_REF] below. In this way, we introduce a cell-like approximation algorithm

µ k (dx) pk (x) dx , with pk (x) = i∈I k μi k λ i k 1 B i k (x) ,
where, for all i ∈ I k , λ i k is the Lebesgue measure of the cell B i k , and μi k is an approximation of the conditional probability

µ i k = P (X t k ∈ B i k | Z k )
. This kind of approximation was proposed by James-Le Gland [START_REF] James | Numerical approximation for nonlinear filtering and finite-time observers[END_REF], for the approximation of nonlinear filters and observers.

Here, the problem is to define, at each time t k , the cells {B i k , i ∈ I k } and the approximate conditional probabilities {μ i k , i ∈ I k }. Among many possible choices, we will focus on the following :

B i k = ξ t k-1 ,t k (B i k-1 ) and μi k = c k R i k μi k-1 , (7)
for all i ∈ I k , where, by definition, I k = I does not depend on k, R i k is the generalized likelihood function for the estimation, given the observation z k , of the parameter i ∈ I such that {X t k ∈ B i k }, that is :

R i k = max x∈B i k Ψ k (x) , (8)
and c k is a normalization constant.

The main points of this section are the following :

(i ) From the numerical point of view, we obtain an estimate of the error

µ k (dx) -pk (x) (dx) ,
in terms of the discretization parameter H, see Theorem 4.1 below.

(ii ) From the point of view of estimating the cell containing the initial condition X 0 ∈ IR m , we show the consistency of the Bayesian parameter in the small noise asymptotics (for a fixed discretization parameter H). This example is also illustrated by the problem already studied at Section 3 for the particle approximation.

The long time asymptotics is more difficult to handle. Actually, in the purely theoretical setup (without approximation) we have recently obtained some results concerning the convergence of the filter to the true value, when the identifiability hypothesis is fulfilled, see Cérou [START_REF] Érou | Long time asymptotics for some dynamical noise free non-linear filtering problems[END_REF][START_REF] Érou | Long time asymptotics for some dynamical noise free non-linear filtering problems, new cases[END_REF].

Problem setting

We consider a particular nonlinear filtering problem where there is no noise input in the state equation. As a result, the only unknown quantity is the initial state of the system.

We can tackle this problem in two ways. Because the only unknown parameter is the initial condition X 0 , we can expect that, in this case, the nonlinear filtering problem reduces to the problem of estimating the parameter X 0 . We can consider the maximum likelihood estimator or the Bayesian estimator, and study the consistency properties, the rate of convergence, etc.

On the other hand, we can study the consequences of this particular setting on the nonlinear filtering equations and their numerical solution.

In the sequel, we study the following model

Ẋt = b(X t ) , X 0 unknown, z k = h(X t k ) + v k , (9)
where

t 0 < t 1 < • • • < t k < • •
• is a strictly increasing sequence of observation times, and {v k , k ≥ 0} is a i.i.d. sequence of centered Gaussian random variables with covariance matrix R. Throughout this paper, we assume for simplicity that the sequence of observation times is uniform, i.e. t k+1 -t k = ∆ for all k = 0, 1, • • •, and we make the following Hypothesis 2.1 The covariance matrix R is non singular. In the case where X 0 is a random variable, we suppose that it is independent of {v k , k ≥ 0}.

Nonlinear filtering approach

Let Z k denotes the σ-algebra generated by the observations up to time t k

Z k = σ(z 1 , • • • , z k ) ,
and suppose that X 0 is a r.v. with probability distribution :

µ 0 (dx) = P (X 0 ∈ dx) .
Our goal is to compute the conditional probability distribution µ k (dx

) = P (X t k ∈ dx | Z k ) of X t k given Z k .
Definition 2.2 (Conditional probability distributions) We introduce the following notation :

(i ) µ k (dx) = P (X t k ∈ dx | Z k ), is the conditional probability distribution of X t k given Z k . (ii ) µ - k (dx) = P (X t k ∈ dx | Z k-1 ), is the conditional probability distribution of X t k given Z k-1 . (iii ) For t k ≤ t ≤ t k+1 , µ k t (dx) = P (X t ∈ dx | Z k ) is the conditional probability distribution of X t given Z k . We have µ k t k = µ k and µ k t k+1 = µ - k+1 .
Proposition 2.3 (Optimal nonlinear filter) The sequence {µ k , k ≥ 0} satisfies a recurrence equation, and the iteration µ k → µ k+1 splits in two steps : prediction, and correction.

Prediction step : From t k to t k+1 , µ k t (dx) satisfies, in a weak sense, the Fokker-Planck equation

∂µ k t ∂t = L * µ k t , (10)
where

L = m i=1 b i ∂ ∂x i
is the partial differential operator associated with the state equation in model [START_REF] Kushner | Numerical methods for stochastic control problems in continuous time[END_REF].

Correction step : At time t k+1 , the a priori information µ - k+1 (dx), is combined with the new observation z k+1 , according to the Bayes formula

µ k+1 (dx) = c k+1 Ψ k+1 (x) µ - k+1 (dx) , (11)
where by definition Ψ k+1 (x) is the likelihood function for the estimation of then parameter X t k+1 ∈ IR m given the observation z k+1

Ψ k+1 (x) = exp -1 2 z k+1 -h(x) 2 R -1 (12)
and c k+1 is a normalization constant.

Proof For t ≥ t k , and for any test function ϕ defined on IR m , we have

ϕ(X t ) = ϕ(X t k ) + t t k Lϕ(X s ) ds ,
and so

E[ϕ(X t )|Z k ] = E[ϕ(X t k )|Z k ] + t t k E[Lϕ(X s )|Z k ] ds , or µ k t , ϕ = µ k , ϕ + t t k µ k s , Lϕ ds .
which proves that {µ k t , t k ≤ t ≤ t k+1 } satisfies the Fokker-Planck equation [START_REF] Le Gland | Application du filtrage non linéaire oprtimal à la restitution d'orbite de transfert[END_REF] in a weak sense.

The correction step is obvious. 2

Parametric estimation approach

We can reformulate the state parameter estimation for the partially observed system [START_REF] Kushner | Numerical methods for stochastic control problems in continuous time[END_REF]. Using the flow of diffeomorphisms ξ s,t (•), we get X t k = ξ 0,t k (X 0 ), for all k ≥ 1, and the observation z k reads

z k = h(ξ 0,t k (X 0 )) + v k . ( 13
)
This is a standard statistical model for the estimation of the unknown parameter X 0 . Actually, we have to choose among trajectories {ξ 0,t (x) , t ≥ 0} for different initial conditions x ∈ IR m at time 0.

Maximum likelihood estimate

The initial condition X 0 (or the state X t at a given time t ≥ 0) is considered as a parameter of IR m , without a priori information. The likelihood function for the estimation of the unknown parameter X 0 in the statistical model defined in [START_REF] Raviart | An analysis of particle methods[END_REF] above, given the observations

{z 1 , • • • , z k } is Ξ k (x) = exp -1 2 k l=1 z l -h(ξ 0,t l (x)) 2 R -1 = k l=1 Ψ l (ξ 0,t l (x)) , (14)
where, for all l = 1, • • • , k, the function Ψ l (•) is defined by [START_REF] Mckeague | Estimation for diffusion processes under misspecified models[END_REF]. The maximum likelihood estimator X 0 is given by :

X 0 ∈ Arg max x∈IR m Ξ k (x) .
Actually, the estimator is not only X 0 , but also the trajectory which satisfies the state equation and starts from X 0 at time 0.

Bayesian estimator

In this section, we have an a priori information on the initial state X 0 , represented as a probability distribution µ 0 (dx) on IR m . This a priori information can be translated, through the state equation and the associated flow of diffeomorphisms, into an a priori information on the state X t k at time t k . We can get an explicit expression for the conditional probability distribution µ k (dx), using the flow of diffeomorphisms ξ s,t (•) associated with the state equation : Proposition 2.4 For any Borel set A ⊂ IR m , we have

µ k (A) = c k ξ -1 0,t k (A) Ξ k (x) µ 0 (dx) ,
where the normalization constant c k is given by

c k = IR m Ξ k (x) µ 0 (dx) , and Ξ k (x) is defined by (14).
Proof First, we translate -via ξ 0,t k (•) -the a priori information on the initial condition X 0 into an a priori information on the state X t k at time t k .

Indeed, X t k = ξ 0,t k (X 0 ) and for any test function ϕ defined on IR m

E[ϕ(X t k )] = E[ϕ(ξ 0,t k (X 0 ))] = ϕ(ξ 0,t k (x)) µ 0 (dx) .
This relation defines the probability distribution µ k 0 (dx) = P (X t k ∈ dx) of the state X t k , in the following way : for any test function ϕ defined on IR m

µ k 0 , ϕ = ϕ(x) µ k 0 (dx) = ϕ(ξ 0,t k (x)) µ 0 (dx) .
Actually, µ k 0 (dx) is the image of the probability distribution µ 0 (dx) under the diffeomorphism ξ 0,t k (•). Equivalently, for any Borel set A ⊂ IR m , we have

µ k 0 (A) = µ 0 (ξ -1 0,t k (A)) .
From the Bayes rule, the conditional probability distribution µ k (dx) of the state X t k , given observations Z k , is given -up to a normalization factor -as the product of the a priori probability distribution µ k 0 (dx) and the likelihood function Ξ k (ξ -1 0,t k (x)) for the estimation of the parameter X t k (or the corresponding initial state ξ -1 0,t k (X t k )), that is

µ k (dx) = c k Ξ k (ξ -1 0,t k (x)) µ k 0 (dx) .
Hence, for any test function ϕ defined on IR m

µ k , ϕ = ϕ(x) µ k (dx) = c k ϕ(x) Ξ k (ξ -1 0,t k (x)) µ k 0 (dx) = c k ϕ(ξ 0,t k (x)) Ξ k (x) µ 0 (dx) ,
and for any Borel set A ⊂ IR m , we have the desired formula. 2

A posteriori probability distribution

The computation of the a posteriori (i.e. given the observations) conditional probability distribution µ k (dx) shows off two interesting particular cases, depending on the form of the probability distribution µ 0 (dx) of the random variable X 0 :

(i ) If µ 0 (dx) has a density with respect to the Lebesgue measure on IR m , then the conditional probability distribution µ k (dx) has also a density.

(ii ) On the other hand, if µ 0 (dx) is a discrete probability distribution (i.e. a linear and convex combination of Dirac measures), then the conditional probability distribution µ k (dx) is also discrete.

Definition 2.5

We define the operator Q k which relates the conditional probability distribution µ k (dx) with the probability distribution µ 0 (dx) of the random variable X 0 :

Q k µ 0 = µ k , that is Q k µ 0 , ϕ = c k ϕ(ξ 0,t k (x)) Ξ k (x) µ 0 (dx) , (15)
for any test function ϕ defined on IR m , or equivalently

Q k µ 0 (A) = c k ξ -1 0,t k (A) Ξ k (x) µ 0 (dx) , for any Borel set A ⊂ IR m .
Proposition 2.6 (µ 0 absolutely continuous) Suppose that the initial probability distribution µ 0 (dx) has a density p 0 (x) with respect to the Lebesgue measure, µ 0 (dx) = p 0 (x) dx. Then µ k (dx) has also a density p k (x) : µ k (dx) = p k (x) dx, defined by

p k (x) = c k [J k (ξ -1 0,t k (x))] -1 Ξ k (ξ -1 0,t k (x)) p 0 (ξ -1 0,t k (x))
, where J k denotes the Jacobian determinant associated with the diffeomorphism ξ 0,t k (•).

Proof By definition of the operator Q k , we have :

µ k , ϕ = Q k µ 0 , ϕ = c k ϕ(ξ 0,t k (x)) Ξ k (x) p 0 (x) dx . Taking ϕ = f • ξ -1 0,t k ,
where f is any test function defined on IR m , we get :

µ k , ϕ = c k f (x) Ξ k (x) p 0 (x) dx . ( 16
)
On the other hand

µ k , ϕ = ϕ(x) p k (x) dx = f (ξ -1 0,t k (x)) p k (x) dx .
The change of variable : x = ξ -1 0,t k (x), gives :

µ k , ϕ = f (x ) p k (ξ 0,t k (x )) J k (x ) dx . (17)
From ( 16) and ( 17) we deduce

p k (ξ 0,t k (x)) J k (x) = c k Ξ k (x) p 0 (x) ,
which gives the probability density function p k (x). 2

Proposition 2.7 (µ 0 discrete) Suppose that the initial probability distribution µ 0 (dx) is a linear combination of Dirac measures :

µ 0 = i∈I a i 0 δ x i 0 ,
where {x i 0 , i ∈ I} and {a i 0 , i ∈ I} are respectively the positions and the weights of the particles. Then the conditional probability distribution µ k (dx) is also a linear combination of Dirac measures :

µ k = i∈I a i k δ x i k , with x i k = ξ 0,t k (x i 0 ) and a i k = c k a i 0 Ξ k (x i 0 ) .
Proof By definition of the operator Q k we have :

µ k , ϕ = Q k µ 0 , ϕ = c k i∈I a i 0 ϕ(ξ 0,t k (x i 0 )) Ξ k (x i 0 ) ,
for any test function ϕ defined on IR m , i.e.

µ k , ϕ = i∈I a i k ϕ(x i k ) .

Particle approximation

In this section, we consider the case where the probability distribution µ 0 (dx) has a density p 0 (x) with respect to the Lebesgue measure on IR m . We introduce µ H 0 (dx), an approximation of µ 0 (dx), as a linear and convex combination of Dirac measures, called particles. The probability distribution µ H k (dx) = Q k µ H 0 (dx) is an approximation of the probability distribution µ k (dx) = Q k µ 0 (dx), and we study the associated approximation error. This kind of approximation was proposed by Raviart [START_REF] Raviart | An analysis of particle methods[END_REF] for deterministic, first order PDE's.

Then, we will see that a coarse approximation, with respect to the covariance matrix R of the observation noise, can produce a non-consistent approximated Bayesian estimator.

Choice of the approximation

Let µ H 0 (dx) be the approximation of the initial probability distribution p 0 (x) dx. We have :

p 0 (x) dx = µ 0 (dx) ∼ µ H 0 (dx) , with µ H 0 = i∈I a i 0 δ x i 0 ,
where {a i 0 , i ∈ I} are the weights of the particles, and {x i 0 , i ∈ I} are the positions of the particles.

First, we fix ε > 0. Then there exist a compact set K ⊂ IR m such that

µ 0 (K ) ≥ 1 -ε .
We introduce a covering of K consisting of bounded and convex Borel sets {B i , i ∈ I} with mutually disjoint interiors, e.g. cubes. We set

K = i∈I B i ⊃ K . A fortiori µ 0 (K) ≥ µ 0 (K ) ≥ 1 -ε . (18)
We define the weights in the following way :

a i 0 = µ 0 (B i ) = B i p 0 (x) dx, a i 0 > 0 . (19)
We can always suppose that a i 0 > 0, otherwise we replace K by K \ B i . Then, we define the position the following way :

x i 0 = 1 a i 0 B i x p 0 (x) dx . ( 20
)
Because B i is convex, we have x i 0 ∈ B i . For all i ∈ I, let δ i be the diameter of the bounded subset B i , and let H be the maximum of the diameters {δ i , i ∈ I}. In particular, for all i ∈ I we get sup

x∈B i |x -x i 0 | ≤ H . (21)
This approximation of µ 0 (dx) has already been considered in Florchinger-Le Gland [START_REF] Florchinger | Particle approximation for first order SPDE's[END_REF].

Fabien Campillo, Frédéric Cérou, François Le Gland, Rivo Rakotozafy

If we use µ H 0 (dx) as an approximation of the initial probability distribution p 0 (x) dx, then the conditional probability distribution p k (x) dx is approximated by µ H k (dx) = Q k µ H 0 (dx), which is a linear combination of Dirac measures, and we get :

p k (x) dx = µ k (dx) ∼ µ H k (dx) , with µ H k = i∈I a i 0 Ξ k (x i 0 ) δ x i k ,
where

x i k = ξ 0,t k (x i 0 ).

Algorithm

We can decompose the particle approximation in two steps

Prediction step : for all i ∈ I x i k+1 = Φ ∆ (x i k ) . (22)
i.e.

x i k+1 is the image of x i k by the diffeomorphism Φ ∆ (•) = ξ t k ,t k+1 (•). Correction step : for all i ∈ I a i k+1 = c k+1 Ψ k+1 (x i k+1 ) a i k , (23)
where c k+1 is a normalization constant and

Ψ k+1 (x) = exp -1 2 z k+1 -h(x) 2 R -1 , x ∈ IR d .

Error estimate

Theorem 3.1 Suppose that b(•) and h(•) are bounded, together with their derivatives up to order 2. Then :

E µ k -µ H k -2,1 = E sup f ∈W 2,∞ | µ k , f -µ H k , f | f 2,∞ ≤ 2 ε + H 2 (C + C r ) .
where H is the largest of the diameters {δ i , i ∈ I}, and r is the smallest eigenvalue of the observation noise covariance matrix R.

First we introduce some notations. We consider the following factorization :

Ξ k (x) = exp -1 2 k l=1 z l 2 R -1 Λ k (x) .
for all x ∈ IR m , which defines Λ k (x). Let P † be the probability measure under which {z k , k ≥ 1} is an i.i.d. sequence of centered Gaussian random variables with covariance matrix R, independent of X 0 .

For all x ∈ IR m , we define the probability measure P x equivalent to P † , with Radon-Nikodym derivative :

dP x dP † Z k = Λ k (x) . ( 24
)
Under the probability measure P x :

z k = h(ξ 0,t k (x)) + v x k ,
for all k ≥ 1, where {v x k , k ≥ 1} is an i.i.d. sequence of centered Gaussian random variables, with covariance matrix R, independent of X 0 .

Finally, the probability measure P satisfies

dP dP † Z k = Λ k (x)µ 0 (dx) = µ 0 , Λ k . (25)
Then, for any test function f defined on IR m

µ k , f = Q k µ 0 , f = f (ξ 0,t k (x)) Ξ k (x) p 0 (x) dx Ξ k (x) p 0 (x) dx = f (ξ 0,t k (x)) Λ k (x) p 0 (x) dx Λ k (x) p 0 (x) dx = µ 0 , g k µ 0 , Λ k , with g k (x) = f (ξ 0,t k (x)) Λ k (x) for all x ∈ IR m . Equivalently : µ H k , f = Q k µ H 0 , f = µ H 0 , g k µ H 0 , Λ k ,
We finally notice that :

µ k , f -µ H k , f = µ 0 , g k µ 0 , Λ k - µ H 0 , g k µ H 0 , Λ k = µ 0 , g k -µ H 0 , g k µ 0 , Λ k - µ H 0 , g k µ H 0 , Λ k µ 0 , Λ k -µ H 0 , Λ k µ 0 , Λ k (26) = E k (f ) -µ H k , f E k (1) ,
where, for any test function

f defined on IR m E k (f ) = µ 0 , g k -µ H 0 , g k µ 0 , Λ k , with g k (x) = f (ξ 0,t k (x)) Λ k (x) for all x ∈ IR m . Hence, it is sufficient to estimate E k (f ), which
is the purpose of the following Lemmas 3.3 and 3.4. But first, we make an error on the initial probability distribution, which we have to evaluate :

Lemma 3.2
The following estimate holds :

µ 0 -µ H 0 -2,1 = sup f ∈W 2,∞ | µ 0 , f -µ H 0 , f | f 2,∞ ≤ ε + 1 2 H 2 .
Fabien Campillo, Frédéric Cérou, François Le Gland, Rivo Rakotozafy Proof Let f be a test function defined on IR m . Taylor expansion of f at point x i 0 reads :

f (x) = f (x i 0 ) + (x -x i 0 ) * f (x i 0 ) +(x -x i 0 ) * 1 0 (1 -u) f [ux + (1 -u)x i 0 ] du (x -x i 0 ) . Moreover : µ 0 , f = K c f (x) p 0 (x) dx + i∈I B i f (x) p 0 (x) dx ,
and µ H 0 , f = i∈I a i 0 f (x i 0 ) .
The difference µ 0 , f -µ H 0 , f satisfies :

µ 0 , f -µ H 0 , f = K c f (x) p 0 (x) dx + i∈I B i [f (x) -f (x i 0 )] p 0 (x) dx = K c f (x) p 0 (x) dx (27) + i∈I B i (x -x i 0 ) * 1 0 (1 -u) f [ux + (1 -u)x i 0 ] du (x -x i 0 ) p 0 (x) dx ,
since by definition (20)

B i (x -x i 0 ) * f (x i 0 ) p 0 (x) dx = 0 , i ∈ I .
Hence we get :

| µ 0 , f -µ H 0 , f | ≤ f ∞,K c K c p 0 (x) dx + 1 2 i∈I f ∞,B i B i |x -x i 0 | 2 p 0 (x) dx .
Moreover, we have

K c p 0 (x) dx = µ 0 (K c ) ≤ ε , (28)
and

B i |x -x i 0 | 2 p 0 (x) dx ≤ H 2 a i 0 ,
according to (18), ( 19) and (21). Hence :

| µ 0 , f -µ H 0 , f | ≤ ε f ∞,K c + 1 2 H 2 f ∞,K ≤ (ε + 1 2 H 2 ) f 2,∞ . 2 Lemma 3.3
The following inequality holds :

E sup f ∈W 2,∞ |E k (f )| f 2,∞ ≤ ε + 1 2 i∈I sup x∈B i E x sup f ∈W 2,∞ |d k (x)| f 2,∞ δ 2 i a i 0 , with g k (x) = f (ξ 0,t k (x)) Λ k (x) and g k (x) = d k (x) Λ k (x) for all x ∈ IR m .
Proof First, according to (25)

E sup f ∈W 2,∞ |E k (f )| f 2,∞ = E sup f ∈W 2,∞ | µ 0 , g k -µ H 0 , g k | f 2,∞ µ 0 , Λ k = E † sup f ∈W 2,∞ | µ 0 , g k -µ H 0 , g k | f 2,∞ .
Moreover, from (27) we get

µ 0 , g k -µ H 0 , g k = K c g k (x) p 0 (x) dx + i∈I B i (x -x i 0 ) * 1 0 (1 -u) g k [ux + (1 -u)x i 0 ] du (x -x i 0 ) p 0 (x) dx .
The following inequality comes up :

| µ 0 , g k -µ H 0 , g k | ≤ K c |g k (x)| p 0 (x) dx + i∈I B i |(x -x i 0 ) 1 0 (1 -u) g k [ux + (1 -u)x i 0 ] du (x -x i 0 )| p 0 (x) dx ≤ f ∞ K c Λ k (x) p 0 (x) dx + i∈I δ 2 i B i 1 0 (1 -u) |d k [ux + (1 -u)x i 0 ]| Λ k [ux + (1 -u)x i 0 ] du p 0 (x) dx ,
where δ i is the diameter of the subset B i . From this estimate, we deduce

sup f ∈W 2,∞ | µ 0 , g k -µ H 0 , g k | f 2,∞ ≤ K c Λ k (x) p 0 (x) dx + i∈I δ 2 i B i 1 0 (1 -u) sup f ∈W 2,∞ |d k [ux + (1 -u)x i 0 ]| f 2,∞ Λ k [ux + (1 -u)x i 0 ] du p 0 (x) dx , hence E † sup f ∈W 2,∞ | µ 0 , g k -µ H 0 , g k | f 2,∞ ≤ sup x∈K c E † [Λ k (x)] K c p 0 (x) dx + 1 2 i∈I δ 2 i sup x∈B i E † sup f ∈W 2,∞ |d k (x)| f 2,∞ Λ k (x) B i p 0 (x) dx . Notice that, for all x ∈ IR m , E † [Λ k (x)] = 1 and E † sup f ∈W 2,∞ |d k (x)| f 2,∞ Λ k (x) = E x sup f ∈W 2,∞ |d k (x)| f 2,∞ ,
according to (24). So, from ( 19) and (28), we get

E † sup f ∈W 2,∞ | µ 0 , g k -µ H 0 , g k | f 2,∞ ≤ ε + 1 2 i∈I sup x∈B i E x sup f ∈W 2,∞ |d k (x)| f 2,∞ δ 2 i a i 0 ,
which completes the proof of the lemma, according to (29).

E x sup f ∈W 2,∞ |d k (x)| f 2,∞ ≤ C + C r where g k (x) = f (ξ 0,t k (x))Λ k (x), and g k (x) = d k (x)Λ k (x)
for all x ∈ IR m , and r is the smallest eigenvalue of the observation noise covariance matrix R.

Proof An explicit computation of the second derivative of g k (x), gives the following expression for the function d k (x) :

d k (x) = f (ξ 0,t k (x)) (ξ 0,t k (x)) 2 + f (ξ 0,t k (x)) ξ 0,t k (x) +2f (ξ 0,t k (x)) ξ 0,t k (x) (log Λ k ) (x) + f (ξ 0,t k (x)) (log Λ k ) (x) +f (ξ 0,t k (x)) [(log Λ k ) (x)] 2 ,
where we suppose for simplicity that m = 1.

Since the functions f (•) and b(•) are bounded together with their derivatives up to order 2, there exists C > 0 such that :

sup f ∈W 2,∞ |d k (x)| f 2,∞ ≤ C[1 + |(log Λ k ) (x)| 2 + |(log Λ k ) (x)|] . Moreover log Λ k (x) = k l=1 z * l R -1 h(ξ 0,t l (x)) -1 2 k l=1 h(ξ 0,t l (x)) 2 R -1 , (log Λ k ) (x) = k l=1 [z l -h(ξ 0,t l (x))] * R -1 h (ξ 0,t l (x)) ξ 0,t l (x) = k l=1 [R -1/2 v x l ] * R -1/2 h (ξ 0,t l (x)) ξ 0,t l (x) , and
(log Λ k ) (x) = k l=1 [z l -h(ξ 0,t l (x))] * R -1 [h (ξ 0,t l (x)) (ξ 0,t l (x)) 2 +h (ξ 0,t l (x)) ξ 0,t l (x)] - k l=1 [h (ξ 0,t l (x)) ξ 0,t l (x)] * R -1 h (ξ 0,t l (x)) ξ 0,t l (x) = k l=1 [R -1/2 v x l ] * R -1/2 [h (ξ 0,t l (x)) (ξ 0,t l (x)) 2 + h (ξ 0,t l (x)) ξ 0,t l (x)] - k l=1 [h (ξ 0,t l (x)) ξ 0,t l (x)] * R -1 h (ξ 0,t l (x)) ξ 0,t l (x)
So we get :

E x |(log Λ k ) (x)| 2 ≤ C r and E x |(log Λ k ) (x)| ≤ C + C r which leads to E x sup f ∈W 2,∞ |d k (x)| f 2,∞ ≤ C + C r . 2
Proof of Theorem 3.1 According to (26), we have

sup f ∈W 2,∞ | µ k , f -µ H k , f | f 2,∞ ≤ sup f ∈W 2,∞ |E k (f )| f 2,∞ + |E k (1)| ≤ 2 sup f ∈W 2,∞ |E k (f )| f 2,∞ .
From estimates proved in Lemmas 3.3 and 3.4, we get

E sup f ∈W 2,∞ | µ k , f -µ H k , f | f 2,∞ ≤ 2ε + H 2 (C + C r ) ,
which proves the theorem. 2

Consistency

The above error estimate shows that it is not sufficient for the discretization step H to be small. It is also necessary for H 2 to be small compared with the smallest eigenvalue r of the observation noise covariance matrix R. The aim of this section is to prove that, if this is not the case, the particle approximation can produce a coarse, but also non-consistent, estimator. To get a better view of this situation, we consider the case where R = rI and r tends to 0. When r > 0, the likelihood function Ξ k (x) for the estimation of the initial value

X 0 satisfies -r log Ξ k (x) = 1 2 k l=1 z l -h(ξ 0,t l (x)) 2 .
The maximum likelihood estimator X 0 is given by

X 0 ∈ Arg max x∈IR m Ξ k (x) .
When r ↓ 0, we get the following limiting expression (Kullback-Leibler information) :

-r log Ξ k (x) -→ K(x, x 0 ) = 1 2 k l=1 h(ξ 0,t l (x 0 )) -h(ξ 0,t l (x)) 2 ,
where x 0 ∈ IR m denotes the true value of the initial condition. We introduce the set

M (x 0 ) = Arg min x∈IR m K(x, x 0 ) = {x ∈ IR m : h(ξ 0,t l (x)) = h(ξ 0,t l (x 0 )) , for all l = 1, • • • , k} M (x 0
) is the set of initial values which, in the limiting deterministic system, cannot be distinguished from the true value x 0 . Obviously, x 0 ∈ M (x 0 ), but the system may be not identifiable, so that M (x 0 ) = {x 0 }. An example of such a system is presented in Lévine-Marino [START_REF] Évine | Constant-speed target tracking via bearings-only measurements[END_REF], where a target with constant speed is tracked with angle measurements only. In this example, the set M (x 0 ) is a one dimensional submanifold.

We have the following consistency result for the maximum likelihood estimator :

d( X 0 , M (x 0 )) → 0 , with probability one, as r ↓ 0 .

The particle approximation described above consists in restricting the parameter set to a finite set G H = {x i 0 , i ∈ I} ⊂ IR m of possible initial values, and nothing can insure that the true value x 0 belongs to G H : this is a mis-specified statistical model, see for example McKeague [START_REF] Mckeague | Estimation for diffusion processes under misspecified models[END_REF]. The maximum likelihood estimator X H 0 is given by

X H 0 ∈ Arg max x∈G H Ξ k (x) .
We define :

M H (x 0 ) = Arg min x∈G H K(x, x 0 ) .
and we get the following consistency result for the maximum likelihood estimator :

d( X H 0 , M H (x 0
)) → 0 , with probability one, as r ↓ 0 .

Usually x 0 ∈ M H (x 0 ) except if x 0 ∈ G H . It can happen that d(x 0 , M H (x 0 )
) is large, in particular when the system is not identifiable. This phenomenon is illustrated in Figure 1, where the set M H (x 0 ) reduce to the single point X H 0 = x 0 , while the set M (x 0 ) corresponds to the continuous curve crossing the true value x 0 .

Cell approximation

For the particle method presented above, we approximate the initial probability distribution µ 0 (dx) by a linear and convex combination of Dirac measures, located at points {x i 0 , i ∈ I}, and we choose among a finite number of possible trajectories starting from these initial conditions.

For the cell approximation, which was first introduced in James-Le Gland [START_REF] James | Numerical approximation for nonlinear filtering and finite-time observers[END_REF], at each time t k , we consider a family of bounded Borel sets {B i k , i ∈ I k } called cells, with mutually disjoint interiors, and we choose between a finite number of composite hypotheses {H i , i ∈ I k }, where for all i ∈ I k , the composite hypothesis H i is {X t k ∈ B i k }. Throughout this section, we assume that the flow of diffeomorphism {ξ s,t (•) , 0 ≤ s ≤ t} associated with ( 9) is explicitly known. In particular, the diffeomorphism

Φ ∆ (•) = ξ t k ,t k+1 (•) is explicitly known.
This assumption may appear quite restrictive : however, there exist interesting problems where this assumption if satisfied. This is the case in the target motion analysis problem [START_REF] Évine | Constant-speed target tracking via bearings-only measurements[END_REF] already mentioned. Numerical results for this problem are presented in Section 6.2.

Choice of the approximation

We introduce the following notation : for all i ∈ I k

µ i k-1/2 = P (X t k ∈ B i k |Z k-1 ) µ i k = P (X t k ∈ B i k |Z k ) .
First, we suppose that the partitions are given at each time, and we present an approximation μk = {μ i k , i ∈ I k } of the discrete probability distribution µ k = {µ i k , i ∈ I k }. We shall later consider the problem of the choice of the partitions.

The computation is done in two steps : prediction and correction.

Prediction

From the discrete probability distribution μk = {μ i k , i ∈ I k } we consider pk (x) dx, the approximation of µ k (dx), given by pk

(x) = i∈I k μi k λ i k 1 B i k (x) , (29)
where λ i k is the Lebesgue measure of the set B i k . Because Φ ∆ (•) is a diffeomorphism, we get : see Figure 2, hence :

{X t k+1 ∈ B i k+1 } = {X t k ∈ Φ -1 ∆ (B i k+1 )} ,
µ i k+1/2 = P (X t k+1 ∈ B i k+1 |Z k ) = P (X t k ∈ Φ -1 ∆ (B i k+1 )|Z k ) = Φ -1 ∆ (B i k+1 ) µ k (dx) .
Using the approximation pk (x) dx defined by (29), we have

µ i k+1/2 j∈I k μj k λ j k Φ -1 ∆ (B i k+1 ) 1 B j k (x) dx = j∈I k μj k λ[B j k ∩ Φ -1 ∆ (B i k+1 )] λ[B j k ]
,

and we define the discrete probability distribution μk+1/2 = {μ i k+1/2 , i ∈ I k+1 } as follows :

μi k+1/2 = j∈I k μj k λ[B j k ∩ Φ -1 ∆ (B i k+1 )] λ[B j k ] . (30)

Correction

From the discrete probability distribution μk+1/2 = {μ i k+1/2 , i ∈ I k+1 } defined by (30), an approximation of µ k+1/2 (dx) = µ - k+1 (dx) is given by pk+1/2 (x) dx, with

pk+1/2 (x) = i∈I k+1 μi k+1/2 λ i k+1 1 B i k+1 (x) , (31)
where λ i k+1 is the Lebesgue measure of the cell B i k+1 . The Bayes formula leads to :

µ i k+1 = P (X t k+1 ∈ B i k+1 |Z k+1 ) = c k+1 B i k+1 Ψ k+1 (x)µ k+1/2 (dx) , INRIA
where c k+1 is a normalization constant.

Using the approximation pk+1/2 (x) dx defined by (31), we have

µ i k+1 c k+1 μi k+1/2 1 λ i k+1 B i k+1 Ψ k+1 (x) dx c k+1 μi k+1/2 R i k+1 ,
where R i k+1 is the generalized likelihood ratio corresponding to the composite hypothesis {X t k+1 ∈ B i k+1 }, that is :

R i k+1 = max x∈B i k+1 Ψ k+1 (x) .
Then, we define the discrete probability distribution μk+1 = {μ i k+1 , i ∈ I k+1 } by :

μi k+1 = c k+1 R i k+1 μi k+1/2 , (32)
where c k+1 is a normalization constant. Let us suppose that the generalized likelihood ratio R i k+1 is explicitly given. This is the case for the target motion analysis problem presented in Section 6.2.

Choice of the partitions

Now we consider the problem of the choice of the partitions B k = {B i k , i ∈ I k } for all k = 0, 1, • • •. As we did above, we fix ε > 0. There exists a compact set K ⊂ IR m such that

µ 0 (K ) ≥ 1 -ε .
We introduce a covering of K consisting of bounded and convex Borel sets {B i 0 , i ∈ I 0 } with mutually disjoint interiors. We set K = i∈I 0

B i 0 ⊃ K . A fortiori µ 0 (K) ≥ µ 0 (K ) ≥ 1 -ε .
We choose I k = I 0 = I for all k = 0, 1, • • •, and for all i ∈ I we set B i k+1 = Φ ∆ (B i k ). This choice leads to :

B j k ∩ Φ -1 ∆ (B i k+1 ) = ∅ if j = i B i k if j = i , so that μi k+1/2 = μi k for all i ∈ I.

Algorithm

The cell approximation consists in computing, for all i ∈ I

μi k+1 = c k+1 R i k+1 μi k , with R i k+1 = exp -1 2 min x∈B i k+1 z k+1 -h(x) 2 R -1
, where the cell B i k+1 = Φ ∆ (B i k ) is the image of the cell B k by the diffeomorphism Φ ∆ (•), and c k+1 is a normalization constant.

Error estimate

Under smoothness assumptions for the coefficients of system (9) we have the following error estimate.

Theorem 4.1 Suppose that b(•) and h(•) are bounded, together with there derivatives up to order 1. Then we get :

E µ k -pk 1 ≤ A k ,
where, for all p ≥ 1, the sequence {A k , k ≥ 0} satisfies the following estimate

A k+1/2 ≤ A k + C H k (33) A k+1 ≤ C p A k+1/2 + C p 1 √ r H β k+1 , (34)
with A 0 ≤ C H 0 . The constant β < 1-d/p could be chosen arbitrarily, H k is the larger of the diameters {δ i k , i ∈ I}, and r is the smallest eigenvalue of the observation noise covariance matrix R.

Proof To prove this result, our goal is to establish, by induction over the index k, the following estimate :

| µ k , f -pk , f | ≤ i∈I sup x∈B i k |f (x)| α i k , (35)
for any test function f defined on IR m . Suppose that the estimate (35) is true : then, for any bounded test function f , we get :

| µ k , f -pk , f | ≤ f ∞ i∈I α i k ,
so :

E µ k -pk 1 = E sup f ∈L ∞ | µ k , f -pk , f | f ∞ ≤ i∈I E[α i k ] = A k ,
which proves the desired error estimate. First, we prove (35) for the initial condition (k = 0). To establish the induction hypothesis, we will study successively the prediction step (from k to k + 1/2), and the correction step (from k + 1/2 to k + 1).

Initial condition Suppose that the probability distribution µ 0 (dx) is absolutely continuous, i.e. µ 0 (dx) = p 0 (x) dx, and has a compact support

K ⊂ i∈I B i 0 . By definition p0 (x) = i∈I μi 0 λ i 0 1 B i 0 (x) with μi 0 = µ 0 (B i 0 ) for all i ∈ I, so that p0 , f = i∈I μi 0 λ i 0 B i 0 f (x) dx .
Hence :

µ 0 , f -p0 , f = i∈I B i 0 f (x) [p 0 (x) - 1 λ i 0 B i 0 p 0 (x ) dx ] dx ,
which leads to the following estimate

| µ 0 , f -p0 , f | ≤ i∈I sup x∈B i 0 |f (x)| 1 λ i 0 B i 0 B i 0 |p 0 (x) -p 0 (x )| dx dx .
The induction hypothesis (35) is proved, with

α i 0 = 1 λ i 0 B i 0 B i 0 |p 0 (x) -p 0 (x )| dx dx ≤ C δ i 0 |p 0 | 1,B i 0 ,
where δ i 0 is the diameter of the cell B i 0 . Notice that :

A 0 = i∈I α i 0 ≤ C H 0 |p 0 | 1 ,
where H 0 denotes the largest of the diameters {δ i 0 , i ∈ I}.

Prediction Introduce the following decomposition (i ) First we have

µ k+1/2 , f -pk+1/2 , f = µ k+1/2 , f -qk+1/2 , f + qk+1/2 , f -pk+1/2 , f ,
µ k+1/2 , f -qk+1/2 , f = µ k , g -pk , g , with g(x) = f (Φ ∆ (x)) for all x ∈ IR m .
From the induction hypothesis (35), we get :

| µ k+1/2 , f -qk+1/2 , f | ≤ i∈I sup x∈B i k |g(x)| α i k .
Notice that : sup

x∈B i k |g(x)| = sup x∈B i k |f (Φ ∆ (x))| = sup x∈B i k+1 |f (x)| , taking into account the choice of B i k = Φ -1 ∆ (B i k+1
). From this we obtain :

| µ k+1/2 , f -qk+1/2 , f | ≤ i∈I sup x∈B i k+1 |f (x)| α i k .
Fabien Campillo, Frédéric Cérou, François Le Gland, Rivo Rakotozafy

(ii ) Moreover : qk+1/2 , f -pk+1/2 , f = i∈I μi k λ i k+1 B i k f (Φ ∆ (x)) [ λ i k+1 λ i k -J ∆ (x)] dx ,
where J ∆ denotes the Jacobian determinant of the diffeomorphism Φ ∆ (•). Hence :

| qk+1/2 , f | -pk+1/2 , f | ≤ i∈I sup x∈B i k+1 |f (x)| μi k λ i k+1 B i k 1 λ i k B i k |J ∆ (x ) -J ∆ (x)| dx dx ≤ i∈I sup x∈B i k+1 |f (x)| μi k C i k ,
with :

C i k = 1 λ i k+1 B i k 1 λ i k B i k |J ∆ (x ) -J ∆ (x)| dx dx ≤ C δ i k .
(iii ) Collecting the above estimates, leads to :

| µ k+1/2 , f -pk+1/2 , f | ≤ i∈I sup x∈B i k+1 |f (x)| [α i k + μi k C i k ] .
Hence, the estimate ( 35) is proved, with :

α i k+1/2 = α i k + μi k C i k .
Notice that :

A k+1/2 = i∈I E[α i k+1/2 ] = A k + i∈I E[μ i k ] C i k ≤ A k + C H k ,
where H k denotes the largest of the diameters {δ i k , i ∈ I}, which proves the estimate (33).

Correction First, we introduce some notations. Consider the following factorization :

Ψ k+1 (x) = exp -1 2 z k+1 2 R -1 Ψ k+1 (x) , (36)
for all x ∈ IR m , which defines Ψ k+1 (x). Let P † k+1 be the probability measure under which

z l = h(X t l ) + v l for all l = 1, • • • , k, and {v 1 , • • • , v k , z k+1 } are i.i.d. centered
Gaussian random variables with covariance matrix R, independent of X 0 . For all x ∈ IR m , we define the probability measure P x k+1 equivalent to P † k+1 , with Radon-Nikodym derivative

dP x k+1 dP † k+1 Z k+1 = Ψ k+1 (x) . ( 37
)
Under the probability measure

P x k+1 z l = h(X t l ) + v l ,
for all l = 1, • • • , k, and

z k+1 = h(x) + v x k+1 , where {v 1 , • • • , v k , v x
k+1 } are i.i.d. centered Gaussian random variables with covariance matrix R, independent of X 0 .

Finally, the probability measure P satisfies

dP dP † k+1 Z k+1 = Ψ k+1 (x)µ k+1/2 (dx) = µ k+1/2 , Ψ k+1 . ( 38
)
From ( 11) and (36), we get

µ k+1 , f = µ k+1/2 , Ψ k+1 f µ k+1/2 , Ψ k+1 = µ k+1/2 , Ψ k+1 f µ k+1/2 , Ψ k+1 ,
whereas, from ( 32) and (36) we get

pk+1 , f = i∈I μi k R i k+1 λ i k+1 B i k+1 f (x) dx i∈I μi k R i k+1 (39)
with : R i k+1 = sup

x∈B i k+1 Ψ k+1 (x) .
We have the equivalent of formula (26), that is :

µ k+1 , f -pk+1 , f = E k+1 (f ) -pk+1 , f E k+1 (1) , (40)
where, for any test function f defined on IR m

E k+1 (f ) = µ k+1/2 , Ψ k+1 f - i∈I μi k R i k+1 λ i k+1 B i k+1 f (x) dx µ k+1/2 , Ψ k+1 .
Then we introduce the following decomposition :

E k+1 (f ) = µ k+1/2 , Ψ k+1 f -pk+1/2 , Ψ k+1 f µ k+1/2 , Ψ k+1 + i∈I μi k λ i k+1 B i k+1 [ Ψ k+1 (x) -R i k+1 ] f (x) dx µ k+1/2 , Ψ k+1 = E k+1 (f ) + E k+1 (f ) ,
and it is sufficient to estimate E k+1 (f ) and E k+1 (f ) separately.

(i ) On one hand :

E k+1 (f ) = µ k+1/2 , g -pk+1/2 , g
where, for all x ∈ IR m

g(x) = Ψ k+1 (x) f (x) µ k+1/2 , Ψ k+1 .
From the induction hypothesis, we get :

|E k+1 (f )| ≤ i∈I sup x∈B i k+1 |g(x)| α i k+1/2 ≤ i∈I sup x∈B i k+1 |f (x)| R i k+1 µ k+1/2 , Ψ k+1 α i k+1/2 .
(ii ) On the other hand :

|E k+1 (f )| ≤ i∈I μi k λ i k+1 B i k+1 [ R i k+1 -Ψ k+1 (x)] |f (x)| dx µ k+1 , Ψ k+1 ≤ i∈I sup x∈B i k+1 |f (x)| μi k λ i k+1 B i k+1 [ R i k+1 -Ψ k+1 (x)] dx µ k+1 , Ψ k+1
(iii ) Collecting the estimates obtained for E k+1 (f ) and E k+1 (f ), leads to :

|E k+1 (f )| ≤ i∈I sup x∈B i k+1 |f (x)| ω i k+1 with ω i k+1 = R i k+1 µ k+1/2 , Ψ k+1 α i k+1/2 + 1 λ i k+1 B i k+1 [ R i k+1 -Ψ k+1 (x)] dx µ k+1/2 , Ψ k+1 μi k . (41)
In particular, we get :

|E k+1 (1)| ≤ i∈I ω i k+1 = ωk+1 . (42)
(iv ) From definition (39), we get

| pk+1 , f | ≤ i∈I sup x∈B i k+1 |f (x)| μi k R i k+1 j∈I μj k R j k+1 . Hence : | pk+1 , f | |E k+1 (1)| ≤ i∈I sup x∈B i k+1 |f (x)| ω i k+1
with :

ω i k+1 = μi k R i k+1 j∈I μj k R j k+1 ωk+1
Notice that :

i∈I ω i k+1 = ωk+1 . (43)
(v ) From (40) and above estimates, we finally get :

| µ k+1 , f -pk+1 , f | ≤ |E k+1 (f )| + | pk+1 , f | |E k+1 (1)| ≤ i∈I sup x∈B i k+1 |f (x)| [ω i k+1 + ω i k+1 ] .
Hence, estimate (35) is proved, with 42) and (43), we have :

α i k+1 = ω i k+1 + ω i k+1 . From (
A k+1 = i∈I E[α i k+1 ] = 2 E[ω k+1 ] . (44)
Therefore, it is sufficient to estimate E[ω i k+1 ] for all i ∈ I, which is the purpose of the three following lemmas. Lemma 4.2 For all i ∈ I, we get

E[ω i k+1 ] = E † k+1 [ R i k+1 ] E[α i k+1/2 ] + E † k+1 [ R i k+1 ] -1 E[μ i k ] . (45)
Proof First, we have :

E[ R i k+1 µ k+1/2 , Ψ k+1 α i k+1/2 ] = E † k+1 [ R i k+1 α i k+1/2 ] = E † k+1 [ R i k+1 ] E[α i k+1/2 ] ,
according to (38), and the independence property of R i k+1 and Z k under the probability measure P † k+1 . With the same arguments, and taking into account the fact that

E † k+1 [ Ψ k+1 (x)] = 1 for all x ∈ IR m , we get E[ 1 λ i k+1 B i k+1 [ R i k+1 -Ψ k+1 (x)] dx µ k+1/2 , Ψ k+1 μi k ] = E † k+1 [ 1 λ i k+1 B i k+1 [ R i k+1 -Ψ k+1 (x)] dx μi k ] = E † k+1 [ 1 λ i k+1 B i k+1 [ R i k+1 -Ψ k+1 (x)] dx] E[μ i k ] = E † k+1 [ R i k+1 ] -1 E[μ i k ] .
Then, the lemma results from definition (41). 2

According to the previous lemma, we just need to estimate

E † k+1 [ R i k+1 ]. Lemma 4.3 For all i ∈ I, we get E † k+1 [ R i k+1 ] ≤ (1 + u p ) p , where u p = E † k+1 [∆ p p ]
1/p and ∆ p = sup

x∈B i k+1 | Ψ 1/p k+1 (x) -Ψ 1/p k+1 (x)| , for any x in B i k+1 .
Proof Let p ≥ 1 be arbitrarily given. Notice that :

[ R i k+1 ] 1/p = sup x∈B i k+1 Ψ 1/p k+1 (x) ≤ Ψ 1/p k+1 (x) + sup x∈B i k+1 | Ψ 1/p k+1 (x) -Ψ 1/p k+1 (x)| ≤ Ψ 1/p k+1 (x) + ∆ p ,
where the definition of ∆ p is given in the statement of the lemma, and x is any point in B i k+1 . From this we deduce :

R i k+1 ≤ ( Ψ 1/p k+1 (x) + ∆ p ) p = p n=0 C n p Ψ n/p k+1 (x) ∆ p-n p .
From the Hölder inequality

E † k+1 [ R i k+1 ] ≤ p n=0 C n p E † k+1 [ Ψ n/p k+1 (x) ∆ p-n p ] ≤ p n=0 C n p E † k+1 [ Ψ rn/p k+1 (x)] 1/r E † k+1 [∆ r (p-n) p ] 1/r . We choose r = p n ≥ 1, so that rn p = 1 and r (p -n) = r r -1 (p -n) = p.
From the fact that E † k+1 [ Ψ k+1 (x)] = 1, we get :

E † k+1 [ R i k+1 ] ≤ p n=0 C n p E † k+1 [∆ p p ] (p-n)/p = p n=0 C n p u p-n p = (1 + u p ) p ,
where

u p = E † k+1 [∆ p p ] 1/p
, which is the desired result. 2

According to the previous lemma, we just have now to estimate u p .

Lemma 4.4 For all i ∈ I, and p ≥ 1, we get

u p ≤ C p 1 √ r [δ i k+1 ] β ,
where β < 1 -d/p, δ i k+1 is the diameter of the cell B i k+1 and r is the smallest eigenvalue of the observation noise covariance matrix R.

Proof For x , x ∈ B i k+1 , Taylor expansion of the function Ψ 1/p k+1 reads :

Ψ 1/p k+1 (x) -Ψ 1/p k+1 (x ) = (x -x ) * 1 0 1 p Ψ 1/p-1 k+1 Ψ k+1 [ux + (1 -u)x ] du = 1 p (x -x ) * 1 0 Ψ k+1 Ψ k+1 Ψ 1/p k+1 [ux + (1 -u)x ] du ,
and we get the following inequality

| Ψ 1/p k+1 (x) -Ψ 1/p k+1 (x )| p ≤ 1 p p |x -x | p 1 0 Ψ k+1 Ψ k+1 p Ψ k+1 [ux + (1 -u)x ] du .
Taking expectation w.r.t. the probability measure P † k+1 , we get from (37)

E † k+1 Ψ 1/p k+1 (x) -Ψ 1/p k+1 (x ) p ≤ 1 p p |x -x | p 1 0 E ux+(1-u)x k+1 Ψ k+1 Ψ k+1 p [ux + (1 -u)x ] du ≤ 1 p p |x -x | p sup x ∈IR m E x k+1 Ψ k+1 Ψ k+1 p (x ) .
Moreover, we have :

Ψ k+1 Ψ k+1 (x) = [z k+1 -h(x)] * R -1 h (x) = [R -1/2 v x k+1 ] * R -1/2 h (x) .
Since the function h(•) and its first derivative are bounded, we get the inequality :

Ψ k+1 Ψ k+1 (x) p ≤ C 1 √ r p |R -1/2 v x k+1 | p ,
where r is the smallest eigenvalue of the covariance matrix R. Taking expectation w.r.t. the probability measure P x k+1 , leads to :

E x k+1 Ψ k+1 Ψ k+1 (x) p ≤ C E |ξ| √ r p ,
where ξ is a centered Gaussian random variable, with identity covariance matrix. From this we derive the following inequality :

E † k+1 | Ψ 1/p k+1 (x) -Ψ 1/p k+1 (x )| p ≤ C p |x -x | √ r p ,
where C p is a constant which depends only on p.

From the Kolmogorov criterion, see for example Kunita [8, Theorem 1.4.1], we get the following inequality :

E † k+1 [∆ p p ] ≤ C p 1 √ r p [δ i k+1 ] βp , with β < 1 -d/p, hence u p ≤ C p 1 √ r [δ i k+1 ] β ,
where δ i k+1 is the diameter of the cell B i k+1 . 2

Fabien Campillo, Frédéric Cérou, François Le Gland, Rivo Rakotozafy

End of the proof of Theorem 4.1 Let us suppose that 1

√ r [δ i k+1 ] β < C. Then E † k+1 [ R i k+1 ] ≤ 1 2 C p , (46)
and

E † k+1 [ R i k+1 ] -1 ≤ 1 2 C p 1 √ r [δ i k+1 ] β . ( 47
)
On the contrary, if the diameter δ i k+1 of the cell B i k+1 is too large, so that the condition 1 √ r [δ i k+1 ] β < C is not fulfilled, we just have to subdivide this cell into subcells with sufficiently small diameters.

In order to prove the estimate (34), and conclude the proof of Theorem 4.1, we just have to combine (44), (45) and the above estimates ( 46) and (47). 2

Consistency

Error estimates seem to indicate that the cell approximation is less accurate than the particle approximation. However, we will show that the cell approximation leads to a consistent estimator of the cell containing the initial condition, as the noise observation covariance matrix tend to zero.

For simplicity, we consider the case where R = rI with r → 0. For r > 0, the generalized likelihood function {Ξ i k , i ∈ I} for the estimation of the cell containing the initial condition X 0 satisfies :

-r log Ξ i k = -r log k l=1 R i l = 1 2 k l=1 min x∈B i 0 z l -h(ξ 0,t l (x)) 2 .
Proposition 4.5 The maximum likelihood estimator î0 is given by :

î0 ∈ Arg min i∈I Ξ i k .
We have the following consistency result for the maximum likelihood estimator : d(î 0 , I(x 0 )) → 0 , with probability one, as r ↓ 0 , where

I(x 0 ) = Arg min i∈I K i (x 0 ) = i ∈ I : min x∈B i 0 h(ξ 0,t l (x 0 )) -h(ξ 0,t l (x)) 2 = 0 , ∀l ≤ k .
I(x 0 ) is the set of initial cells which contain, for each observation time t l , at least one initial condition which, in the limiting deterministic system, cannot be distinguished from the true value x 0 , based on the observation available at time t l only. Obviously, i 0 ∈ I(x 0 ) if x 0 ∈ B i 0 0 .

Proof When r ↓ 0, we get the following limiting expression :

-r log Ξ i k -→ K i (x 0 ) = 1 2 k l=1 min x∈B i 0 h(ξ 0,t l (x 0 )) -h(ξ 0,t l (x)) 2 ,
where x 0 ∈ IR m denotes the true value of the initial condition. 2

INRIA

Numerical implementation

At this level it is necessary to define what we exactly want to get as an output of the (approximate) nonlinear filter. This question is closely related to the question of graphical output. The conditional mean and covariance are quite poor estimators (i.e. functions of the conditional probability distribution). In fact, the concept of confidence region is much more meaningful. This concept is defined below for an absolutely continuous probability distribution. Therefore, in the case of the particle approximation, where we have only a discrete approximation of the conditional probability distribution, it will be necessary, in a first step, to get an approximate conditional density.

Here we suppose that I k = I 0 = I for all k. Given the conditional density p k we compute an approximation of Dα k as follows. Let

Confidence regions

p max = max x∈IR n p k (x)
We divide the interval [0, p max ] into K sub-intervals [η i , η i+1 ], with η 0 = 0 and η K = p max , and we set

D i k = {x ∈ IR n : p k (x) ≥ η i }.
Then by putting

i 0 = max{i : µ k (D i k ) ≥ α} ,
we obtain B i 0 k as an approximate confidence region.

Particle approximation

The nonlinear filter

Numerically, formula (23) is not well suited. Indeed, it may happen (especially if the covariance of the observation noise is small) that the likelihood Ψ k+1 (x i k+1 ) is very small for some i ∈ I, hence numerically zero, and then we will have a l = 0 for all l ≥ k + 1. A way to overcome this difficulty is to use a logarithmic transformation, i.e. we set

l i k = log a i k . (48)

Density reconstruction

In order to get graphical outputs with confidence regions it is necessary to transform the conditional probability distribution given as a linear combination of Dirac measures, into a conditional density. We propose the following simple algorithm.

Initialization

First, we take a regular bounded subdomain D of IR n and we usually take as initial condition particles uniformly distributed on D : we define G n δ as

G n δ = D ∩ IR n δ
where IR n δ is the n-dimensional grid with mesh size δ. We put the particles on the nodes of the grid, and we take uniform initial weights, i.e.

a i 0 = 1 N , for all i = 1, • • • , N , where N = Card G n δ .
This initialization is simple and does not assume any a priori information.

Construction of an approximate conditional density

In order to compute confidence regions, we need to transform the discrete conditional probability distribution into an absolutely continuous probability distribution. We take a partition of the domain D : {A j , j ∈ J}. Then by putting

λ j k = 1 µ(A j ) i∈I a i k 1 A j (x i k ) , ∀j ∈ J, ∀k ∈ IN .
we obtain an approximation of the conditional density by pk

(x) = j∈J λ j k 1 A j (x)
for all x ∈ D. For simplicity we choose a regular partition of the domain D, see Figure 3 for a simple 2-dimensional example.

Cell approximation

For the particle method presented above, we approximate the initial probability distribution µ 0 (dx) by a linear and convex combination of Dirac measures, located at points {x i 0 , i ∈ I}, and we choose among a finite number of possible trajectories starting from these initial conditions.

For the cell approximation, at each instant t k , we consider a family of bounded Borel sets {B i k , i ∈ I} called cells, with mutually disjoint interiors, and we choose between a finite number of composite hypotheses {H i , i ∈ I}, where for all i ∈ I, the composite hypothesis Prediction

H i is {X t k ∈ B i k }. For a given family of partitions B k = {B i k , i ∈ I}, we compute the following approxima- tions : μi k-1/2 P (X t k ∈ B i k |Z k-1 ) μi k P (X t k ∈ B i k |Z k ) for all i ∈ I,
B i k+1 = Φ ∆ (B i k ) , i ∈ I . (50)
Correction

μi k+1 = c k+1 R i k+1 μi k+1/2 ,
where c k+1 is a normalization constant and

R i k+1 = max x∈B i k+1 Ψ k+1 (x) . Algorithm Let νi k = log μi k , i ∈ I Initialization : μi 0 ← µ 0 (B i 0 ) , i ∈ I . then νi k+1 ← ρ i k+1 + νi k , i ∈ I l * ← max i∈I νi k+1 , νi k+1 ← νi k+1 -l * , i ∈ I μi k+1 ← exp{ν i k+1 } , i ∈ I c k+1 ← i∈I μi k+1 μi k+1 ← 1 c k+1 μi k+1 , i ∈ I with ρ i k+1 = -1 2 min x∈B i k+1 z k+1 -h(x) 2 R -1 = log R i k+1 . (51)
where B i k+1 = Φ ∆ (B i k ), see Section 4.1 for the choice of the partitions, and c k+1 is a normalization constant.

The reason for introducing the normalization l * is to prevent νi k+1 from taking increasingly large negative values.

In fact, this algorithm is not completely practical, because in general :

• The prediction step (50) is not explicit.

• The computation of the maximum in (51) is also not explicit.

To get a practical algorithm, we have to consider real case situations. In the first example, see Section 6.1 below, the flow is an isometry, so the two previous points are explicit. In the second example, see Section 6.2 below, an additional approximation is introduced in the prediction step so as to get a partition into parallelepipeds. As a consequence, we have to use the general formula (30) instead. The minimization on the parallelepipeds is then explicit.

Parallel computing

The algorithms presented in the previous sections are well adapted to parallel computation.

In the case of a noisy state equation, it is standard to use upwind finite difference schemes, so as to insure the positivity of the solution, and to obtain a probabilistic interpretation for the approximation, see Kushner-Dupuis [START_REF] Kushner | Numerical methods for stochastic control problems in continuous time[END_REF]. For a multi-dimensional system, this method requires many points of discretization, and is quite slow on a sequential computer.

The same problem arises for the particle / cell methods. We can think of implementing these methods on a vector supercomputer, but one can easily see that our algorithms are much more adapted to parallel than vector processing.

As the main part of the algorithms is local on each discretization point, or involves only neighbors, it can be computed on all nodes at the same time in a very efficient way.

6 Numerical results

Example 1 : particle vs. cell approximation

Let us consider the following test problem :

Ẋt = A X t , X 0 unknown, z k = h(X t k ) + v k , ( 52
)
where X is a process which take values in IR 2 and {v k , k ≥ 1} is a sequence of i.i.d. N (0, σ 2 ) random variables.

x 1

x 2

In this problem, the state of the system has a circular motion. Here, the matrix A is given by :

A = 0 -1 1 0
Let ∆ = t k+1 -t k denotes the time between two observations, then the flow associated with the system (52) is :

Φ ∆ = cos ∆ -sin ∆ sin ∆ cos ∆
In this case, we suppose that we have access only to information concerning the distance between the mobile and the origin (0, 0). The observation function is given by :

h(X t k ) = x 2 1 (t k ) + x 2 2 (t k ) , X t k = x 1 (t k ) x 2 (t k ) .
The numerical values used for simulation are :

• Initial probability distribution : uniform law over [-5, 5] × [-5, 5],

• Time step between two measurements : ∆ = 0.1.

• Standard deviation of the observation noise : σ = 0.5.

• We take the following values for the initial conditions :

X 0 = 4 0 and X = 2 √ 2 2 √
2 .

• Rectangular and uniform grid on [-5, 5] × [-5, 5].

• Number of particles : 20 × 20 and 100 × 100. In this section we deal with the following problem : we want to estimate motion parameters of a target in a plane (typically the surface of the sea). The only informations we have about the target comes from bearing measurements made from a moving observer. We suppose that the target has a constant velocity, and that we have discrete time measurements, see Figure 14. Then the state equation is

Ẋ1 t = X 3 t Ẋ2 t = X 4 t Ẋ3 t = 0 Ẋ4 t = 0
where X 1 t and X 2 t denote the coordinates of the target in the plane, and X 3 t , X 4 t denote the two components of the (constant) target velocity vector. The observations are given by

z k = arctan X 1 t k -X o,1 t k X 2 t k -X o,2 t k + v k
where X o,1 t and X o,2 t denote the coordinates of the observation platform at time t, and {v k , k ≥ 0} is a sequence of i.i.d. N (0, r 2 ) random variables. Since this system has a noisefree state equation, it will be solved by the particle method.

Numerical results

The parameters used for the simulation are the following :

• t max = 3600, t k = k∆t, ∆t = 12, times of measurements, • v p = 4, speed of the observation platform,

• α 1 = -45, α 2 = 135, α 3 = -45, directions of the observation platform, taken at times i × t max 3 , i = 0, 1, 2,

• σ = 1, standard deviation of the observation noise.

For the four-dimensional discretization, we have:

• x min = -3000, x max = 5000,

• y min = 20000, y max = 50000,

• v min • the number of particles in each direction is 32.

Remark 6.1 The use of the Connection Machine, with the C * programming language, has shown the interest of parallel programming for nonlinear filtering and optimal control. When we began to implement our algorithms it was one of the most efficient parallel computer available, and C * a good high level language on this particular machine. Now, with the new generation of parallel computers, for example the CM-5, and new programming languages, like CM-Fortran, which is a primary implementation of the future High Performance Fortran (HPF), we hope to solve this type of problems faster and for higher dimensions.

Remark 6.2 (Cell approximation) For the cell approximation approach presented in Section 5.3, we must explicit the prediction step (50). In this example we choose to introduce an additional approximation so as to get, after the prediction step, a partition into parallelepipeds.

Prediction step

In this example, we make the following choice for the partition : we start with a regular initial partition {B i 0 ; i ∈ I} of a bounded subset of IR 4 . We suppose that for each i ∈ I, the cell B i 0 is a parallelepiped of center x i 0 , with constant side length δ p along the p-th coordinate, p = 1, • • • , 4. Then at step k + 1 we make the following choice : let x i k+1 = Φ ∆ (x i k ) denotes the center of the parallelepiped B i k+1 with side length δ p along the p-th coordinate, p = 1, • • • , 4. We use the formula (30) :

μi k+1/2 = j∈I μj k λ[B j k ∩ Φ -1 ∆ (B i k+1 )] λ[B j k ]
(λ Lebesgue measure) .

In this case the computation of λ[B j k ∩ Φ -1 ∆ (B i k+1 )] is explicit.

Lemma 3 . 4

 34 If f (•), b(•) and h(•) are bounded, together with their derivatives up to order 2, then there exist C > 0 and C > 0 such that :

Figure 1 :

 1 Figure 1: Example of non-consistency of the Bayesian estimator given by the particle approximation.

Figure 2 :

 2 Figure 2: Cell approximation.

where qk+1/ 2

 2 (x) dx is the image of the probability distribution pk (x) dx under the diffeomorphism Φ ∆ (•), i.e. qk+1/2 , f = f (Φ ∆ (x)) pk (x) dx , for any test function f defined on IR m .

First, we define

  the concept of confidence region of level α ∈ [0, 1]. We denote by µ k the probability distribution with density p k w.r.t. the Lebesgue measure on IR n . Definition 5.1 A confidence region of level α ∈ [0, 1] is a domain Dα k of IR n , with µ kprobability α and least Lebesgue measure, that is Dα k ∈ Arg min ⊂ IR n ; µ k (D) ≥ α} .

APPROXIMATIONFigure 3 :

 3 Figure 3: Construction of an approximate density from a particle set.

Figure 4 :

 4 Figure 4: Conditional density function, at time t = 10, cell approximation, 100 × 100 points grid, observation noise variance 0.01.

Figure 5 :

 5 Figure 5: Conditional density function, at time t = 10, particle approximation, 100 × 100 points grid, observation noise variance 0.01.

Figure 6 :

 6 Figure 6: Cell approximation, at time t = 10, 100 × 100 points grid, observation noise variance 0.01.

Figure 7 :

 7 Figure 7: Cell approximation, at time t = 10, 100 × 100 points grid, observation noise variance 0.1.

Figure 8 :

 8 Figure 8: Cell approximation, at time t = 10, 20 × 20 points grid, observation noise variance 0.1.

Figure 9 :

 9 Figure 9: Cell approximation, at time t = 10, 20 × 20 points grid, observation noise variance 1.

Figure 10 :

 10 Figure 10: Particle approximation, at time t = 10, 100 × 100 points grid, observation noise variance 0.01.

Figure 11 :

 11 Figure 11: Particle approximation, at time t = 10, 100 × 100 points grid, observation noise variance 0.1.

Figure 12 :

 12 Figure 12: Particle approximation, at time t = 10, 20 × 20 points grid, observation noise variance 0.1.

Figure 13 : 1 . 6 . 2

 13162 Figure 13: Particle approximation, at time t = 10, 20 × 20 points grid, observation noise variance 1.

•Figure 14 :

 14 Figure 14: Target motion analysis with bearings-only measurements

Figure 15 :

 15 Figure 15: Position (left) and velocity (right) marginals at the beginning of filtering (t=480s)

RR n • 2567

Fabien Campillo, Frédéric Cérou, François Le Gland, Rivo Rakotozafy INRIA

Fabien Campillo, Frédéric Cérou, François Le Gland, Rivo Rakotozafy

Files

data input file simulation input file ini-density input file filter output file

Reading data

Read N , t 0 in file data Read x i , a i , i = 1, . . . , N in file ini-density

, N in file filter

End do

Table 1: Particle approximation algorithm -the nonlinear filter. Formula (23) reduces to

where

The algorithm is presented in Table 1.

Minimization (51) in the correction step For each i ∈ I, let {P i, , 0 ≤ ≤ 16} denote the corners of the four-dimensional paralleleliped B i k+1 , and let (x 1 i, , x 2 i, , x 3 i, , x 4 i, ) denote the components of the -th corner. The observation depends only on the first two components. So, let

Then

In fact, in formulas (53), minimization and maximization are done on 4 points because only the first two components on the corners are concerned. Remark 6.3 Performance analysis of the particle approximation applied to target motion analysis has been studied in [START_REF] Bernier | Parallélisme de données et filtrage non linéaire -analyse de performance[END_REF]. In this paper we compare the performances obtained on different kind of computers (standard computers, super computers and parallels computers).