Silvère Bonnabel

Jean-Jacques Slotine

Particle observers for state estimation and adaptation in deterministic systems with random piecewise constant inputs

In this paper we develop a methodology for state estimation of partially contracting systems inspired by the particle filter (PF). If a system is partially contracting, in the sense that once some state variables are assumed known and viewed as inputs, the remaining reduced system is contracting, then only those variables need to be estimated. Indeed, once those variables are known, the remaining variables are asymptotically recovered owing to the contraction properties. To estimate the state of such a system, we thus suggest to use a PF approach to estimate only the non-contracting part, whereas the other part is naturally recovered. The methodology is applied to a nontrivial tracking problem. The trajectories of objects one seeks to track, such as aircraft and marine vehicles, typically consist of smooth sections with large, but infrequent, unpredictable changes. As a result, their motion is well modeled by nonlinear deterministic ordinary differential equations driven by piecewise constant random control inputs. Following our methodology, we combine a type of PF called the variable rate particle filter (VRPF) with a bank of observers assumed to have convergence properties. Using the notion of virtual system, we devise a hybrid variable rate particle observer, that uses the particles to detect changes, and the observers for state estimation in between changes.

I. INTRODUCTION

Particle filters have attracted a lot of attention over the past two decades, with the advent of powerful computers, see e.g., [START_REF] Doucet | On sequential Monte Carlo sampling methods for bayesian filtering[END_REF]. Nowadays, they are able to execute at many hundreds of Hertz for realistic problems, as object tracking (with early implementations dating back to [START_REF] Moral | Particle resolution and non-linear signal processing with radar/sonar application[END_REF]), automotive state estimation, and parameter learning. Virtually all practical implementations of the particle filter resort to the Rao-Blackwellized particle filter (RBPF): the state consists of two groups variables (r, z), such that given r the distribution of variable z may be estimated analytically, typically via a Kalman filter. It thus suffices to use particles to estimate r, and, using Bayes' theorem the whole distribution of (r, z) is then estimated. This reduces the number of particles and improves estimation accuracy [START_REF] Doucet | Rao-Blackwellised particle filtering for dynamic bayesian networks[END_REF].

Building upon this idea, we propose a methodology of particle observers (PO) for partially contracting systems, i.e., systems with state (r, z), and such that if r(t) is known and viewed as an input, the reduced system governing z(t) is deterministic and is contracting, see [START_REF] Wang | On partial contraction analysis for coupled nonlinear oscillators[END_REF]. Along the lines of the RBPF methodology, the rationale is to use particles to estimate the statistical distribution of variable r, whereas z may be asymptotically recovered owing to contraction properties [START_REF] Lohmiller | On contraction analysis for nonlinear systems[END_REF]. This leads to a parsimonious representation of the state, and makes the particle filtering problem more tractable and the filter possibly more accurate.

Intuitively, the main difference is as follows. In the RBPF, r(ω,t) and z(ω,t) are stochastic processes, and if r(•,t) is known, then the distribution of z(•,t) is supposed to be analytically computable. In our approach, z(•,t) is deterministically related to r(•,t) owing to the conditionally deterministic assumption. Thus, given r we need to estimate a value for z(t). This is where contraction properties allow us to derive an approximation ẑ(t) of this "true" value z(t) instead of an entire distribution, making the use of an observer more suited to estimating the z part of the state.

We then tackle a tracking application. Trajectories of aircraft and marine vehicles are in fact smooth and well described by continuous time ordinary differential equations (ODE) driven by constant inputs, but subject to abrupt changes at each maneuver. Based on this approach to modeling the tracking task, S. Godsill with various co-authors proposed the variable rate particle filter [START_REF] Godsill | Variable rate particle filters for tracking applications[END_REF], [START_REF] Bunch | Dynamical models for tracking with the variable rate particle filter[END_REF], [START_REF]Particle smoothing algorithms for variable rate models[END_REF], a particle filter well suited to piecewise deterministic models.

On the one hand, it is well known that particle filters have trouble converging in the presence of low process noise, and are subject to divergence when dealing with (even partially) deterministic models, see e.g., [START_REF] Kantas | On particle methods for parameter estimation in state-space models[END_REF]. In particular, to the authors best knowledge, the stability of the variable rate particle filter has never been studied. On the other hand, the field of nonlinear observers, and notably the theory of contraction [START_REF] Lohmiller | On contraction analysis for nonlinear systems[END_REF], provide powerful tools from dynamical systems theory to achieve stable estimation for noise-free models, but are not suited to tasks such as change detection. This prompts an application of our methodology of particle observers for partially contracting systems. We consider deterministic systems with inputs jumping randomly at unknown times, and run many observers in parallel (the particles), each assumed to be contracting in the absence of jumps. This way, we obtain a hybrid filter which is able to detect the jumps, and to reconverge to the true state between jumps. The theory is made rigorous using the notion of virtual systems introduced in [START_REF] Wang | On partial contraction analysis for coupled nonlinear oscillators[END_REF], see also [START_REF] Jouffroy | Methodological remarks on contraction theory[END_REF]. Our approach may be related to the flow particle filter [START_REF] Daum | Nonlinear filters: beyond the Kalman filter[END_REF], and the feedback particle filter [START_REF] Yang | Feedback particle filter[END_REF]. Some preliminary ideas can be found in our recent paper [START_REF] Bonnabel | Particle observers for contracting dynamical systems[END_REF].

Section II is devoted to the particle observer methodology. Section III considers the piecewise deterministic estimation problem, and the variable rate particle filter is introduced. In Section IV, assuming we have at our disposal an observer with some convergence properties, a hybrid strategy that combines it with a particle filter is proposed.

II. PARTICLE FILTERS AND PARTICLE OBSERVERS

A. The Rao-Blackwellized particle filter (RBPF) Consider a discrete Markov process r n having initial distribution p(r 0 ) and transition kernel p(r n | r n-1 ). The variable r n is supposed hidden, and assume we have observations in the form of a random variable y n at time n correlated with r n . Moreover, assume the observations to be conditionally independent given the hidden variables. The goal of discrete time filtering is to estimate the hidden variables at best from the observations, i.e., to ideally compute

p(r n | y 1:n ), where y 1:n = {y 1 , • • • , y n }, (1)
or more generally p(r 1:n | y 1:n ), a task known as smoothing.

Assume now, that observations involve another process z n , such that p(z n | y 1:n , r 1:n ) may be analytically evaluated, typically using a Kalman filter.

If we may approximate the desired distribution using a set of weighted particles, i.e., p(r

1:n | y 1:n ) ≈ ∑ N 1 w (i) n δ r (i) 1:n
(r 1:n ), then the best estimate based on past observations y 1:n of any quantity F(z n , r n ) in the sense of minimal mean squared error may be approximated as:

E p(z n ,r n |y 1:n ) (F(z n , r n )) ≈ N ∑ i=1 w (i) n E p(z n |y 1:n ,r (i) 1:n ) F(z n , r (i) n )
where we used p(z n , r

1:n | y 1:n ) = p(z n | r 1:n , y 1:n )p(r 1:n | y 1:n ).
This gives rise to the celebrated Rao-Blackwellized Particle Filter (RBPF). Algorithm 1 provides its basic version.

Algorithm 1 RBPF with prior sampling Draw N particles from the prior initial distribution p(r 0 )

for i = 1 • • • N do
Sample from the prior

r (i) n ∼ p(r n | r (i)
n-1 ), and let r (i)

1:n = r (i) n , r (i) 1:n-1
Evaluate and update weights

w (i) n = p(y n | y 1:n-1 , r (i) 1:n ) w (i) n-1 Normalize weights w(i) n = w (i) n /[∑ j w ( j) n ] -1 Our estimate for value F(z n , r n ) given past observations is ∑ i w(i) n E p(z n |y 1:n ,r (i) 1:n ) F(z n , r (i) n ) .
Resample if necessary, i.e., multiply/delete particles to obtain N random samples with weights equal to 1/N. end for

B. Particle observers for partially contracting systems

In this paragraph, we summarize the results of [START_REF] Bonnabel | Particle observers for contracting dynamical systems[END_REF]. Consider a noisy dynamical system of the form

d dt z = f (z, r) (2)
d dt r = g(r) + w(t) (3)
where f , g are smooth maps, and w(t) process noise. Definition 1: If equation ( 2) is a contraction when r(t) is considered as a known input, system (2)-( 3) is said to be a partial contraction [START_REF] Wang | On partial contraction analysis for coupled nonlinear oscillators[END_REF]. This implies, see [START_REF] Lohmiller | On contraction analysis for nonlinear systems[END_REF], that if r(t) is known and viewed as an input of (2), then (2) has the property that for any two trajectories z 1 (t), z 2 (t) we have:

d(z 1 (t), z 2 (t)) ≤ d(z 1 (0), z 2 (0))e -2λt , (4)
with λ the contraction rate, and d a well-chosen distance. Assume discrete time noisy measurements are available:

y n = h(z(t n ), r(t n )) +V n at times t 0 < t 1 < t 2 < • • • , (5)
with V n unknown independent identically distributed random variables with density l, i.e., p(y | z, r) = l(yh(z, r)).

Note that, by denoting (r n , z n ) := (r(t n ), z(t n )) one may want to apply the RBPF methodology to dynamics (2)-(3) with observations [START_REF] Lohmiller | On contraction analysis for nonlinear systems[END_REF]. However, the usual assumption underlying the RBPF that p(z n | y 1:n , r 1:n ) may be analytically evaluated does not hold. It may yet be approximated under the partial contraction assumption introduced in Definition 1. Indeed, if we call ẑ(t) a solution of (2) associated to some given trajectory {r(t)} t≥0 , then asymptotically we have p(z(t) | {r(s)} s≤t ) ≈ δ (z(t)ẑ(t)), since (4) implies exponential convergence of ẑ(t) to z(t). Thus, letting (ẑ(t) (i) , r(t) (i) ) be a solution to the stochastic differential equations (SDE) (2)-( 3), we have the approximations:

p(z(t) | y 1:n , {r(s) (i) } s≤t ) ≈ δ (z t -ẑ(t) (i) ), t ≥ t n (6)
E p(z(t)|y 1:n ,{r(s

) (i) } s≤t ) F(z(t), r(t) (i) ) ≈ F(z(t) (i) , r(t) (i) ), p(y n | y 1:n-1 , {r(s) (i) } s≤t n ) ≈ l(y n -h(ẑ(t n ) (i) , r(t n ) (i) )). ( 7
)
An application of the RBPF methodology to (2)-(3) with observations (5) under the partial contraction assumption and the subsequent approximations ( 6)-(7) thus yields the particle observer (PO) Algorithm 2.

Algorithm 2

The PO with prior sampling Draw N particles from the prior initial distribution p(z 0 , r 0 )

for i = 1 • • • N do Sample (z(t n ) (i) , r(t n ) (i)
) by numerically integrating the stochastic differential equations ( 2)-(3) from t n-1 to t n . Evaluate and update weights

w (i) n = l(y n -h(z(t n ) (i) , r(t n ) (i) )) w (i) n-1 .

Normalize weights w(i)

n = w (i) n /[∑ j w ( j)
n ] -1 . Our estimate for value F(z, r) given past observations is

∑ i w(i) n F(r(t) (i) , z(t) (i) ).
In particular the state estimate is ∑ i w(i) n (r(t) (i) , z(t) (i) ).

Resample if necessary. end for

III. STATE ESTIMATION FOR DETERMINISTIC SYSTEMS WITH RANDOM JUMPS

The reader might wonder to what extent it is realistic to consider systems of the form (2)-( 3) where noise affects only [START_REF] Doucet | Rao-Blackwellised particle filtering for dynamic bayesian networks[END_REF]. First, it may serve as a meaningful approximation, if for instance (2) is subject to small noise. Then, and above all, it may readily be applied to various realistic systems. One example consists of nonlinear piecewise deterministic systems subject to abrupt unpredictable jumps. These are relevant to tracking applications, and in the remainder of this paper we apply the ideas of Section II to such systems. Various other perspectives are offered in Section V.

A. Original considered model

In this section we review the conditionally deterministic inference problem addressed by the variable rate particle filter [START_REF] Godsill | Variable rate particle filters for tracking applications[END_REF], see also [START_REF]Particle smoothing algorithms for variable rate models[END_REF] for a somehow more concise presentation. We consider a deterministic dynamical model with continuous-time state x(t) ∈ R n , and parameter (or control input) u ∈ R m . We assume that changepoints occur at random times

{τ 0 = 0, τ 1 , • • • , τ k , • • • }, each with associated parame- ters {u 0 , u 1 , • • • , u k , • • • }. The pairs (τ k , u k ) form elements of a marked point process. Consider some discrete measurement times t 0 ,t 1 , • • • ,t n , • • • , along with associated measurements y 0 , y 1 , • • • , y n , • • • of the form y n = h(x(t n )) +V n (8)
with V n a random variable with known distribution. By letting l denote the probability density function of V n we have

p(y n | x(t n )) = l(y n -h(x(t n ))). ( 9
)
The objective is to estimate from partial and noisy measurements of x(t), the changepoints, the associated parameters, and the state x(t) in real time. Let θ n represent the past and present variables

θ n = {τ 1 , • • • , τ j , u 1 , • • • , u j : 0 ≤ τ j < t n ≤ τ j+1 }.
Let us define the following quantities

K(t) = max{k : τ k < t}, K n = K(t n )
Note that, we could alternatively write K n = K(θ n ) as

K n = max{i : τ i ∈ θ n } (10)
Processes τ k and u k are assumed to be Markov, and the transitions to be given by

p(τ k | τ k-1 ) = p τ (τ k -τ k-1 ), p u (u k | u k-1 )
with p τ a probability density such that p τ (τ) = 0 for τ < 0. The state evolution is assumed to be governed by an ODE which depends on the changepoints and the corresponding parameters as follows:

d dt x(t) = f (x(t),t, u K(t) ). (11)
This means that the state is governed by an ODE driven by parameters that remain constant over random time intervals. Let us denote φ the flow of the ODE ẋ = f (x,t, u), that is,

φ (x 1 , u,t 1 ,t 2 )
is the solution at t 2 to the latter ODE with fixed parameter u and initial condition x 1 at t 1 . We have:

x(t) = φ (x(τ K n ), u K n , τ K n ,t), τ K n < t ≤ τ K n +1 . (12)
By choosing t = τ K n +1 , this equation specifies the state at the next changepoint, and by reiterating to the following changepoints, we get the whole future trajectory. This proves that the knowledge of the process θ n and the initial state x(0) uniquely specifies the state x(t). For simplicity, we will systematically assume x(0) to be known. This is not a very restrictive assumption, as the aim of tracking is to keep track of x(t) despite the random changepoints, and after a few jumps this knowledge about the initial state is forgotten.

Recalling [START_REF] Doucet | On sequential Monte Carlo sampling methods for bayesian filtering[END_REF], the goal of the variable rate particle filter is to compute an approximation to the posterior

p(θ n | y 1:n ) (13)
and x(t) as a byproduct, through [START_REF] Yang | Feedback particle filter[END_REF], i.e. by integrating [START_REF] Daum | Nonlinear filters: beyond the Kalman filter[END_REF].

B. The variable rate particle filter

The bootstrap variable rate particle filter (VRPF) is a particle filter that computes an approximation

p(θ n | y 1:n ) = N ∑ j=1 w ( j) n δ θ ( j) n (θ n ) (14)
to the posterior [START_REF] Bonnabel | Particle observers for contracting dynamical systems[END_REF]. In its simplest bootstrap form, the piecewise deterministic VRPF is implemented along the lines of Algorithm 1, see [START_REF]Particle smoothing algorithms for variable rate models[END_REF], yielding Algorithm 3. 23)-(24) for each particle 1 ≤ j ≤ N, variable χ j (t) tends to move towards true x(t) even if the corresponding tentative couple (τ j , ν j ) is misidentified. The philosophy of particle filtering, that has been compared to Darwin's "survival of the fittest", consists in evaluating weights w j for each particle indicating its fit to measured data. Particles with low weights are discarded, and particles with high weight are multiplied to populate the likely regions of the space, yielding quicker convergence. Here particle i corresponds to a tentative (τ i , ν i ) being close to actual (τ n+1 , u n+1 ) and thus χ i (t) rapidly converges to the true x(t). Particle i is thus likely to be selected and multiplied, whereas particle j having less adequate (τ j , ν j ) might be suppressed.

V. PERSPECTIVES

Building upon the notion of virtual systems [START_REF] Wang | On partial contraction analysis for coupled nonlinear oscillators[END_REF], we proposed a novel method which bridges the gap between nonlinear observer theory and particle filters for a broad class of practical applications.

We anticipate the method might prove useful in various contexts. One such instance is target tracking using one particular model of the nonlinear deterministic models with random piecewise constant inputs proposed in [START_REF] Bunch | Dynamical models for tracking with the variable rate particle filter[END_REF]. For a 3D target model inspired by [START_REF] Bunch | Dynamical models for tracking with the variable rate particle filter[END_REF], an invariant extended Kalman filter based on the work [START_REF] Barrau | The invariant extended Kalman filter as a stable observer[END_REF] was shown in [START_REF] Pilté | Tracking the Frenet-Serret frame associated to a highly maneuvering target in 3D[END_REF] to yield good tracking performance, and might be used as a building block in a particle observer approach.

The method can also be applied to adaptive estimation of deterministic systems parameterized by an unknown vector selected from a finite set of known vectors. The switch is randomly connected to one of the vectors and may either remain there throughout the duration of the process, or change at random times associated to a prior temporal distribution. A typical example application would be Simultaneous Localization and Mapping (SLAM) with a vector accounting for the multi-hypotheses that stem from perceptual ambiguity, as e.g. in [START_REF] Kanitscheider | Training recurrent networks to generate hypotheses about how the brain solves hard navigation problems[END_REF]. To address such questions one could combine the contracting observer of [START_REF] Tan | Analytical SLAM without linearization[END_REF], based on the construction of simple synthetic measurements for SLAM, with the particle approach developed in this article. Another possible application is contact detection in robotics, see also Section IV-D. In the recent paper [START_REF] Bledt | Contact model fusion for event-based locomotion in unstructured terrains[END_REF] for instance, leg state switching is based on a periodic phase-based state schedule, which provides a prior on times of change of a contact parameter s between stance (s = 0) and swing (s = 1).

Disturbance observers [START_REF] Bledt | Contact model fusion for event-based locomotion in unstructured terrains[END_REF] with convergence properties for known s might be run in parallel using our approach.

The methodology of particle observers might also prove useful in the context of differentially positive systems [START_REF] Forni | Differentially positive systems[END_REF], where the state is known to converge to a Perron-Frobenius curve: the use of particles might allow us to identify accurately where the state is located on the curve.

  S. Bonnabel is with Mines ParisTech, PSL Research University, Centre for Robotics, 60 bd Saint-Michel, 75006 Paris, France silvere.bonnabel@mines-paristech.fr J.J. Slotine is with the Department of Mechanical Engineering and the Department of Brain and Cognitive Sciences, Massachusetts Institute of Technology, Cambridge, MA 02139 USA, jjs@mit.edu

Fig. 1 .

 1 Fig. 1. Graphical illustration of Algorithm 4 with x(t) ∈ R 2 . Plain stars indicate changing points at unknown times τ 1 , τ 2 , • • • , each being associated with unknown parameter u i . The particle observer of Algorithm 4 keepsexploring the space of jump times τ and parameter u by constantly creating tentative couples (τ j , ν j ) with associated variable χ j (t) referred to as particles, and makes use of a great number N of such particles. Owing to the assumed convergence properties of observer (23)-(24) for each particle 1 ≤ j ≤ N, variable χ j (t) tends to move towards true x(t) even if the corresponding tentative couple (τ j , ν j ) is misidentified. The philosophy of particle filtering, that has been compared to Darwin's "survival of the fittest", consists in evaluating weights w j for each particle indicating its fit to measured data. Particles with low weights are discarded, and particles with high weight are multiplied to populate the likely regions of the space, yielding quicker convergence. Here particle i corresponds to a tentative (τ i , ν i ) being close to actual (τ n+1 , u n+1 ) and thus χ i (t) rapidly converges to the true x(t). Particle i is thus likely to be selected and multiplied, whereas particle j having less adequate (τ j , ν j ) might be suppressed.

VI. ACKNOWLEDGEMENTS

We thank Prof. Godsill for useful discussions, while the first author was visiting the University of Cambridge as a fellow of Sidney Sussex College. The authors would also like to thank Marion Pilté and Frédéric Barbaresco from the company Thales Air Systems for fruitful discussions.

Algorithm 3 Variable rate particle filter

For all 1 ≤ j ≤ N, initialize changepoint sequence as θ

n-1 ) for each particle j. Using equation [START_REF] Jouffroy | Methodological remarks on contraction theory[END_REF] to get K ( j) n , we calculate the associated x(t) at current measurement time t n :

Evaluate and update weight w

using [START_REF] Kantas | On particle methods for parameter estimation in state-space models[END_REF], and normalize weights to get w( j) n .

Resample if necessary. end for

The seeked distribution is approximated by [START_REF] Barrau | The invariant extended Kalman filter as a stable observer[END_REF]. This yields ∑ j w( j) n x ( j) n as an estimate for the state x(t n ). end for

Let us now comment on some of the steps. Regarding the sampling step, assume that we want to sample τ K n given τ

. This means there has been exactly K ( j) n-1 changepoints before t n-1 . We can sample then the next changepoint time as τ

, then we let τ denote the latter time and we let K

n-1 . θ does not need to be sampled anymore until the first measurement time t k , k > n to become larger than τ. Or τ

< t n , then we integrate the flow [START_REF] Yang | Feedback particle filter[END_REF] to obtain x ( j) (τ

), and reiterate until the changepoint time becomes larger than t n .

Regarding the weight update, particles are sampled using the prior and the weights are classically defined as w

n-1 , see Algorithm 1. As x(0) is known (and thus considered as a fixed parameter), given the sequence θ ( j) n of changepoint times and parameters, the state corresponds to the integration of system [START_REF] Daum | Nonlinear filters: beyond the Kalman filter[END_REF], that is x ( j) n as defined by [START_REF] Pilté | Tracking the Frenet-Serret frame associated to a highly maneuvering target in 3D[END_REF]. Using that x

n ) indeed, as y n is independent of all the other variables conditionally on x n .

C. Possible caveats in a partially deterministic context

It is well known that particle filters are not efficient for systems governed by deterministic dynamics, typically for fixed parameter inference. Actually, they are prone to divergence, see e.g., [START_REF] Kantas | On particle methods for parameter estimation in state-space models[END_REF]. In the context of the VRPF, though, the state space is continually explored thanks to the jumps. But asymptotic convergence may necessitate a large number of particles to get full coverage. Another possible limitation of PF is in the presence of low measurement noise. Since p(y n | x n ) is then peaked, virtually all weights drop to zero after one measurement, leading in turn to estimation failure.

IV. A PARTICLE OBSERVER BASED APPROACH

In the present section, we consider the model of Section III-A, and propose a particle observer, that addresses some of the issues raised in Section III-C.

Up to a few subtleties related to marked point processes in place of discrete-time Markov processes, Algorithm 3 is close in essence to the standard Algorithm 1, since θ = (τ, u) plays the role of hidden variable r, and x plays the role of auxiliary variable z, see notably the discussion after Algorithm 3. As a result, one may want to pursue the methodology of particle observers described at Section II-B for this application. Our idea is that, indeed, if the sequence of jump times τ 1 , τ 2 , • • • were known, then an asymptotic observer might be devised to estimate the remaining variables, that is, x(t) and input u(t). And under contraction properties of the observer's dynamics, approximations akin to ( 6)-( 7) may then be used. To formalize this approach, we will need to split the observations in two. One type of observation will be used by the observer, and the other type by the particle filter to assess the weights. We will assume to have discrete the noisy measurements [START_REF]Particle smoothing algorithms for variable rate models[END_REF], but also some auxiliary continuous-time noise-free measurements of the form:

Note that the two sources of information are assumed independent, and h and h may bear a quite different structure.

For instance, z(t), albeit noise free, may be of very small dimension compared to x(t).

A. Considered problem and further assumptions

To recap, the considered system is

where u K (t) are parameters jumping at random times. Besides, let us consider the following alternative system

which corresponds to the noise free continuous-time system ( 17)-( 19) between two jumps (as the parameters are fixed between jumps). Furthermore, assume we have at our disposal an asymptotic observer for the latter system of the form:

having asymptotic convergence properties for arbitrary initial condition ( x(0), û(0)). Typically, a desirable property would be exponential contraction [START_REF] Lohmiller | On contraction analysis for nonlinear systems[END_REF], leading to

On the one hand, this observer can readily be used on the original problem ( 17)-( 19). This will lead to a solution that is unlikely to diverge (provided average time between consecutive jumps is larger than the convergence time of the observer), but in the transitory phases that follow changepoints, it will yield inaccurate estimates until ( x(t), û(t)) has re-converged to (x(t), u K (t)).

On the other hand, the variable rate particle filter cannot readily account for the information provided by z(t) of eq. ( 19), since noise free measurements (as they appear in the theory of observers), or very small measurement noise, lead to failure of particle filters as explained in Section III-C. Besides, owing to the fact the dynamics are piecewise deterministic, the variable rate particle filter based solely on the information brought by y n of eq. ( 18) will never be able to guarantee asymptotic convergence to the true parameters, even if the time between jumps is large (as the particles u ( j) remain fixed between two jumps). For those reasons, we propose a hybrid strategy, combining both approaches.

B. Proposed hybrid strategy

Our idea is to devise a variable rate particle filter for system (17)-( 18), and to use the side partial information [START_REF] Forni | Differentially positive systems[END_REF] to help the sampled parameters u ( j) (and thus the state) move towards the "true" parameters and state between consecutive changepoints. Crucially, this can be done by resorting to the notion of virtual systems [START_REF] Wang | On partial contraction analysis for coupled nonlinear oscillators[END_REF], [START_REF] Jouffroy | Methodological remarks on contraction theory[END_REF], i.e., by applying the variable rate methodology to an auxiliary virtual system.

Consider indeed the following virtual system associated to the observer

and rewrite those dynamics as

Building upon the latter systems, we let z(t) = h(x(t)) as before, and we define auxiliary dynamics (recall that

that is, we integrate (28) between two consecutive changepoints, and, at each changepoint at time τ k , the variable ν jumps to the corresponding parameter value u k . Letting φ denote the flow of (28), i.e., φ (χ 1 , ν 1 ,t 1 , z (t 1 ,t 2 ) ,t 2 ) is the solution at t 2 to (28) being initialized on (χ 1 , ν 1 ) at t 1 , and where we let z (t 1 ,t 2 ) = ( h(x(s)) t 1 ≤s≤t 2 . Our auxiliary virtual system can be alternatively defined by

(30) that is, we recover an auxiliary system to which the variable rate methodology of Section III-B can be applied. This yields the proposed Algorithm 4.

Algorithm 4 Hybrid variable rate particle observer

For all 1 ≤ j ≤ N, initialize changepoint sequence as θ

n-1 ) for each particle j. Using equation [START_REF] Jouffroy | Methodological remarks on contraction theory[END_REF] to get K ( j) n , we calculate the associated χ(t), ν(t) at current measurement time t n :

Evaluate weight w

n-1 and normalize weights to get w( j) n . Resample if necessary. end for ∑ j w

n is an estimate for x(t n ). end for Using the flow of auxiliary system (26)-( 27) in place of the flow of original system [START_REF] Tan | Analytical SLAM without linearization[END_REF] to estimate the state of the original system through a VRPF, may be justified through the notion of virtual system, see Section IV-C below.

Remark 1: Particle filters have drawn criticism regarding their execution time. In this respect, our method is particularly efficient, owing to the use of observers as building blocks, each being very cheap to implement, as for each particle we need to integrate an ODE and weight computation is also cheap since p(y | χ) = l(yh(χ)), see [START_REF] Kantas | On particle methods for parameter estimation in state-space models[END_REF]. This is in contrast to the standard RBPF where Kalman filters are commonly used to estimate the z variable, see Section II.

C. Discussion and formal results

System ( 17)-( 18)-( 19) is stochastic, since it is driven by the process ω → (τ k (ω), u k (ω)), and ω → V n (ω) is a random variable. Notably, z(t) is a random process ω → z(t, ω). Based on the fact that in the VRPF methodology the initial condition x(0) is assumed to be known, the following result is noteworthy:

Lemma 1: Consider system ( 17)-( 19) on the one hand, and system (29) on the other. Then, if χ(0) = x(0), we have χ(t) = x(t) at all times t ≥ 0.

Proof: It suffices to note that, since z(t) = h(x(t)), we have z(t) -h(x(t)) = 0, and the pair (x(t), u K(t) ) is thus a solution to the stochastic differential equation (29). This property lies at the heart of the use of virtual systems for observer design [START_REF] Wang | On partial contraction analysis for coupled nonlinear oscillators[END_REF], [START_REF] Jouffroy | Methodological remarks on contraction theory[END_REF]. In terms of our particle filter approach, this means Algorithm 4 may be viewed as a (bootstrap variable rate) particle filter for the original problem, given that initial conditions coincide χ(0) = x(0). Indeed, in this case both dynamical systems coincide as explained in the proof of Lemma 1. Normally, Algorithm 4 should use observations related to its own latent variables, i.e., of the form ỹn = h(χ(t n )) + V n , whereas it evaluates the weights based on the measurements y n = h(x(t n )) + V n instead. However, it is legitimate to use y n as measurements for the auxiliary system since we know from Lemma 1 that at all times χ(t) = x(t) and thus ỹn = y n indeed.

However, Lemma 1 does not hold under an arbitrary small error in measurements described by [START_REF] Forni | Differentially positive systems[END_REF], or for χ(0) = x(0). But those are precisely the cases where our observerbased approach applies. Indeed, by assuming observer (23)-(24) to have convergence properties, each simulated particle (χ ( j) (t), ν ( j) (t)) obtained through φ in Algorithm 4 tends to move closer to (x(t), u K(t) ) between jumps. As a result, particles tend to be less dispersed around x(t) and to yield more precise and robust estimates, see Fig. 1 for an illustration.

D. Generalizations

Note that Eq. ( 17) being quite general, our approach can be used to enhance existing adaptive observers to track more rapidly variations in some unknown parameters u subject to abrupt changes.

Another possibility is to let u denote a model index, and then τ denotes the unknown times of switching.

Besides, the assumption (21) of constant u(t) between jumps is not necessary. Between jumps, generic dynamical behavior for u(t) may in fact be assumed as long as it is contracting or may be tracked by an observer. In this case u K(t) in ( 17) shall be replaced with u K(t) (t) and (21) by the ad hoc differential equation d dt u(t) = g(u(t),t).


Mohammad Noormohammadpour

Ajitesh Srivastava

Cauligi S Raghavendra

Ming Hsieh

On Minimizing the Completion Times of Long Flows over Inter-Datacenter WAN

Keywords: Routing, Flow Completion Time, Traffic Engineering, Software Defined Networking, Inter-Datacenter

Long flows contribute huge volumes of traffic over inter-datacenter WAN. The Flow Completion Time (FCT) is a vital network performance metric that affects the running time of distributed applications and the users' quality of experience. Flow routing techniques based on propagation or queuing latency or instantaneous link utilization are insufficient for minimization of the long flows' FCT. We propose a routing approach that uses the remaining sizes and paths of all ongoing flows to minimize the worst-case completion time of incoming flows assuming no knowledge of future flow arrivals. Our approach can be formulated as an NP-Hard graph optimization problem. We propose BWRH, a heuristic to quickly generate an approximate solution. We evaluate BWRH against several real WAN topologies and two different traffic patterns. We see that BWRH provides solutions with an average optimality gap of less than Furthermore, we show that compared to other popular routing heuristics, BWRH reduces the mean and tail FCT by up to 1.46× and 1.53×, respectively.

I. INTRODUCTION

Recently, dedicated inter-datacenter networks have been used by multiple organizations to connect dozens of their datacenters such as Google's B4 [START_REF] Jain | B4: Experience with a globally-deployed software defined wan[END_REF], [START_REF] Hong | B4 and after: Managing hierarchy, partitioning, and asymmetry for availability and scale in google's softwaredefined wan[END_REF], Facebook's Express Backbone [3], and Microsoft's Global WAN [4]. These networks are owned or leased by one organization, are relatively small with tens of edges, and can be managed in a logically centralized manner, for example, using frameworks such as SDN [START_REF] Mckeown | Software-defined networking[END_REF]. This opens up new opportunities for global networkwide optimizations by combining the knowledge of traffic generated at the datacenters and control over traffic forwarding in the network. Therefore, in this paper, we revisit the wellknown flow routing problem over inter-datacenter networks.

We focus on long flows which carry tremendous volumes of data over inter-datacenter networks [START_REF] Jain | B4: Experience with a globally-deployed software defined wan[END_REF], [START_REF] Kandula | Calendaring for wide area networks[END_REF]. They are usually generated as a result of replicating large objects such as search index files, virtual machine migration, and multimedia content. For instance, over Facebook's Express Backbone, about 80% of flows for cache applications take at least 10 seconds to complete [START_REF] Roy | Inside the Social Network's (Datacenter) Network[END_REF]. Besides, the volume of inter-datacenter traffic for replication of content and data, which generates many long flows, has been growing at a fast pace [3].

In general, flows are generated by different applications at unknown times to move data across the datacenters. Therefore, we assume that flows can arrive at the inter-datacenter network at any time and no knowledge of future flow arrivals. Every flow is specified with a source, a destination, an arrival time, and its total volume of data. The Flow Completion Time (FCT) of a flow is the time from its arrival until its completion.

We focus on minimizing the completion times of long flows which is a critical performance metric as it can significantly affect the overall application performance or considerably improve users' quality of experience. For example, in cloud applications such as Hadoop, moving data faster across datacenters can reduce the overall data processing time. As another example, moving popular multimedia content quickly to a regional datacenter via replication allows improved user experience for many local users. To attain this goal, routing and scheduling need to be considered together which can lead to a complex discrete optimization problem. In this paper, we only address the routing problem, that is, choosing a fixed path for an incoming flow given the network topology and the currently ongoing flows while making no assumptions on the traffic scheduling policy. We focus on single path routing which mitigates the undesirable effects of packet reordering.

A variety of metrics have been used for path selection over WAN including static metrics such as hop count and interface bandwidth, and dynamic metrics such as end-to-end latency which is a function of propagation and queuing latency, and current link bandwidth utilization [START_REF] Kodialam | Online multicast routing with bandwidth guarantees: a new approach using multicast network flow[END_REF], [9]. While these metrics are effective for routing of short flows, they are insufficient for improving the completion times of long flows as we will demonstrate. Over inter-datacenter WAN where end-points are managed by the organization that also controls the routing [START_REF] Jain | B4: Experience with a globally-deployed software defined wan[END_REF], [3], [4], one can use routing techniques that differentiate long flows from short flows and use flow properties obtained from applications, including flow size information, to reduce the completion times of long flows.

In this context, our paper makes the following contributions:

• Assuming no knowledge of future flow arrivals and no constraints on the network traffic scheduling policy, we propose to minimize the worst-case completion time of every incoming flow given the network topology, the currently ongoing flows' paths, and their remaining number of data units. For any given scheduling policy, we route the flows to minimize the worst-case flow completion time. We refer to this routing approach as the Best Worstcase Routing (BWR).

• BWR aims to select a minimum weight path for every incoming flow where a path's weight is defined as the total number of remaining data units of all the ongoing flows that have a common edge with the path. It can be shown that BWR is NP-Hard and finding an optimal solution requires examining all existing paths between the two ends of an incoming flow. We develop a heuristic, called BWRH, to quickly compute a route.

• We run extensive simulations to compare BWRH's performance with that of several other routing heuristics, including popular ones. We show that BWRH improves the mean and tail completion times by up to 1.46× and 1.53×, respectively, given various flow size distributions and scheduling policies. We also show that over multiple topologies and with different traffic patterns, BWRH's optimality gap is, on average, below 0.25%.

II. SYSTEM MODEL

We consider a general network topology with bidirectional links and equal capacity of one for all edges and assume an online scenario where flows arrive at unknown times in the future and are assigned a fixed path as they arrive. Each flow is divided into many equal size pieces (e.g., IP datagrams) which we refer to as data units. We also assume knowledge of the flow size (i.e., number of a flow's data units) for the new flow and the remaining flow size for all ongoing flows. Given an index i, every flow F i is defined with a source s i , a destination t i , an arrival time α i , and a total volume of data V i . In addition, each flow is associated with a path P i , a finish time β i which is the time of delivery of its last data unit, and a completion time c i = β i -α i . Finally, at any moment, the total number of remaining data units of F i is V r i ≤ V i . Similar to multiple existing inter-datacenter networks [START_REF] Jain | B4: Experience with a globally-deployed software defined wan[END_REF], [3], [4], we assume the availability of logically centralized control over the network routing. A controller can maintain information on the currently ongoing long flows with their remaining data units and perform routing decisions for an incoming long flow upon arrival.

We employ a slotted timeline model where at each timeslot a single data unit can traverse any path in the network. In other words, we assume a zero propagation and queuing latency which we justify by focusing only on long flows. Given this model, if multiple flows have a shared edge, only one of them can transmit during a timeslot. We say two data units are competing if they belong to flows that share a common edge. Depending on the scheduling policy that is used, these data units may be sent in different orders but never at the same time. Also, if two flows with pending data units use non-overlapping paths, they can transmit their data units at the same time if no other flow with a common edge with either one of these flows is transmitting at the same timeslot.

III. BEST WORST-CASE ROUTING

We aim to reduce long flows' completion times with no assumption on the scheduling policy for transmission of data units. To achieve this goal, we propose the following routing technique referred to as Best Worst-case Routing (BWR):

Problem 1. Given a network topology G(V, E) and the set of ongoing flows F F F = {F i , 1 ≤ i ≤ N }, we want to assign a path P N +1 to the new flow F N +1 so that the worst-case completion time of F N +1 , i.e., max(c N +1 ) is minimized. Assuming no knowledge of future flows and given the described network model, since only a single data unit can get through any edge per timeslot, the worst-case completion time of a flow happens when the data units of all the flows that share at least one edge with the new flow's path go sequentially and before the last data unit of the new flow is transmitted. Therefore, Problem 1 can be reduced to the following graph optimization problem which aims to minimize the number of competing data units with F N +1 .

Problem 2. Given a network topology G(V, E) where every edge e ∈ E is associated with a set of flows F F F e (that is, e ∈ P i , ∀F i ∈ F F F e ), the set of ongoing flows F F F = {F i , 1 ≤ i ≤ N }, and an incoming flow F N +1 , we want to find a minimum weight path P N +1 where the weight of any path P from s N +1 to t N +1 is computed as follows:

W P = {1≤i ≤ N | F i ∈ {∪ e∈ P F F F e } } V r i (1)
Proposition 1. Assuming no knowledge of future flow arrivals, P N +1 selected by solving Problem 2 minimizes the worst-case completion time of F N +1 regardless of the scheduling policy used for transmission of data units.

Proof. P N +1 is chosen to minimize the maximum number of data units ahead of F N +1 given the knowledge of ongoing flows' remaining data units which minimizes the worst-case β N +1 , that is, the maximum number of timeslots the last data unit of F N +1 has to wait before it can be sent. Since α N +1 is fixed, this minimizes max(c N +1 ).

Example: Consider the scenario shown in Figure 1. A new flow F 4 with 3 data units has arrived and has two options of sharing an edge with F 1 that has 4 remaining data units (path 1) or sharing edges with {F 2 , F 3 } which have a total of 6 remaining data units (path 2). Our approach tries to minimize the worst-case completion time of F 4 given ongoing flows. If path 1 is chosen, the worst case completion time of F 4 will be 7 while with path 2 it will be 9 and therefore, the logically centralized network controller will select path 1 for F 4 . The worst-case completion times are not affected by the scheduling policy and are independent of it. Also, the fact that 2 has three common edges with path 2 and F 3 has two common edges with path 2 does not affect the worst-case completion time of F 4 on path 2.

IV. BWR HEURISTIC (BWRH)

The path weight assignment used in Problem 2 is not edgedecomposable. Finding a minimum weight path for

F N +1 is Algorithm 1: BWRH Input: F N +1 , G(V, E), P i , V r i , 1 ≤ i ≤ N Output: P N +1
K ← #hops on the minimum hop path from s N +1 to t N +1 ; W K min ← Weight of the minimum weight path from s N +1 to t N +1 with at most K hops by examining all such paths; repeat NP-Hard and requires examining all paths from s N +1 to t N +1 . 1We propose a fast heuristic here, called BWRH, that finds an approximate solution to Problem 2. Algorithm 1 shows our proposed approach to finding a path P N +1 for F N +1 . At every iteration, the algorithm finds the minimum weight path from s N +1 to t N +1 with at most K hops by computing the weight of every such path according to Eq. 1. The algorithm starts by searching all the minimum hop paths from s N +1 to t N +1 and finding the weight of the minimum weight path among such paths. It then increases the number of maximum hops allowed (i.e., K) by one, extending the search space to more paths. This process continues until the weight of the minimum weight path with at most K hops is the same as K -1, i.e., there is no gain while increasing the number of hops.

K ← K + 1; Compute W K min ; until W K min ≥ W K-
The termination condition used in BWRH may prevent us from searching long paths. Therefore, if the optimal path is considerably longer than the minimum hop path, it is possible that the algorithm terminates before it reaches the optimal path. Let us call the optimal path P o and the path selected with our heuristic P h . The optimality gap, defined as

W P h -W Po W Po
, is highly dependant on the number of remaining data units of ongoing flows. We find that the worst-case optimality gap can be generally unbounded. However, it is highly unlikely, in general, for the optimal path to be long as having more edges increases the likelihood of sharing edges with more ongoing flows which increases the weight of the path. We will later confirm this intuition through empirical evaluations and show that BWRH provides solutions with an average optimality gap of less than a quarter of percent.

V. APPLICATION TO REAL NETWORK SCENARIOS

We discuss how BWRH can be used to find a path for an incoming flow on a real network assuming a uniform link capacity. We can use the same topology as the actual topology as input to BWRH. Since we focus on long flows for which the transmission time is significantly larger than both propagation and queuing latency along existing paths, it is reasonable to ignore their effect in routing (hence the assumption that these values are zero in §II). Next, assuming that all data units are of the same size, we can use the total number of remaining bytes per ongoing flow in place of the number of remaining data units as it does not affect the selected path. In practice, some data units may be smaller than the underlying network's MTU, which for the long flows with many data units, has minimal effect on the selected path. Once BWRH selects a path, the network's forwarding state is updated accordingly to route the new flow's traffic, for example, using SDN [START_REF] Jain | B4: Experience with a globally-deployed software defined wan[END_REF], [START_REF] Kandula | Calendaring for wide area networks[END_REF].

In general, network traffic is a mix of short and long flows. Since our proposal targets the long flows, routing of short flows will not be affected and could be done considering the propagation and queuing latency. Incoming long flows can be routed according to the knowledge of current long flows while ignoring the effect of short flows. Improving the completion times of long flows in a network with a mix of short and long flows is part of the future work.

VI. EVALUATIONS

We considered two flow size distributions of light-tailed (Exponential) and heavy-tailed (Pareto) and considered Poisson flow arrivals with the rate of λ. We also assumed an average flow size of µ data units with a maximum of 500 data units along with a minimum size of 2 data units for the heavy-tailed distribution. We considered the scheduling policies of First Come First Serve (FCFS), Shortest Remaining Processing Time (SRPT) and Fair Sharing based on max-min fairness [START_REF] Bertsekas | Data networks[END_REF].

Topologies: GScale [START_REF] Jain | B4: Experience with a globally-deployed software defined wan[END_REF] with 12 nodes and 19 edges, AGIS [11] with 25 nodes and 30 edges, ANS [12] with 18 nodes and 25 edges, and Cogent [13] with 197 nodes and 243 edges. We assumed bidirectional edges with a uniform capacity of 1 data unit per time unit for all of these topologies.

Schemes: We considered three schemes besides BWRH. The Shortest Path (Min-Hop) approach simply selects a fixed shortest hop path from the source to destination per flow. The Min-Max Utilization approach selects a path that has the minimum value of maximum utilization across all paths going from the source to the destination. This approach has been extensively used in the traffic engineering literature [START_REF] Kandula | Calendaring for wide area networks[END_REF], [START_REF] Kodialam | Online multicast routing with bandwidth guarantees: a new approach using multicast network flow[END_REF]. The Shortest Path (Random-Uniform) selects a path randomly with equal probability across all existing paths which are at most one hop longer than the shortest hop path.

BWRH's Optimality Gap: In Figure 2 we compute the optimality gap of solutions found by BWRH over three different topologies and under two traffic patterns. The optimal solution was computed by taking into account all existing paths and finding the minimum weight path on topologies of GScale, AGIS, and ANS. We also implemented a custom branch and bound approach which would require less computation time with a small number of ongoing flows (i.e., < 20 in our setting) and an intractable amount of time for a large number of ongoing flows (i.e., > 30 in our setting). According to the results, the average gap is less than 0.25% over all experiments. We could not perform this experiment on larger topologies as computing the optimal solution would take an intractable amount of time.

Effect of Scheduling Policies:

In Figure 3 we fixed the flow arrival rate to 1 and mean flow size to 50 and tried various scheduling policies under the Cogent topology which is much larger than GScale, ANS, and AGIS. All simulations were repeated 10 times and the standard deviation for each instance has been reported. The minimum value normalizes each group of bars. We see that BWRH is consistently better than other schemes regardless of the scheduling policy used. We can also see that compared to each other, the performance of other schemes varies considerably with the scheduling policy applied. To quantify, BWRH provides between 1.18× to 1.53× better tail completion times than the other schemes across all scenarios on average. We also observe up to 1.46× better mean completion times compared to other schemes across all scheduling policies on average.

Running Time: We implemented Algorithm 1 in Java using the JGraphT library. To exhaustively find all paths with at most K hops, we used the class AllDirectedPaths in JGraphT. We performed simulations while varying λ from 1 to 10 and µ from 5 to 50 over 1000 flow arrivals per experiment which covers both lightly and heavily loaded regimes. We also experimented with all the four topologies pointed to earlier, both traffic patterns of light-tailed and heavy-tailed, and all three scheduling policies of FCFS, SRPT, and Fair Sharing. The maximum running time of Algorithm 1 was 222.24 milliseconds, and the average of maximum running time across all experiments was 27 milliseconds. This latency can be considered negligible given the time needed to complete long flows once they are routed.

VII. CONCLUSIONS AND FUTURE WORK

We proposed a new technique for routing based on flow size information to reduce flow completion times. Accordingly, the online routing problem turns into finding a minimum weight path on the topology from the source to the destination where the weight is computed by summing up the number of remaining data units of all the flows that have a common edge with the path. Since this is a hard problem, we proposed a fast heuristic with small average optimality gap. We then discussed how information from a real network scenario could be used as input to the proposed network model to find a path on the actual network for an incoming flow. In the future, we would like to study networks with non-uniform link capacity, multipath routing, and the effect of inaccurate flow size information on routing performance.

Fig. 1 .

 1 Fig. 1. Example of routing a new flow F 4

Fig. 2 .

 2 Fig. 2. BWRH's optimality gap for λ = 10 and µ = 50 computed for 1000 flow arrivals.

Fig. 3 .

 3 Fig. 3. Online routing techniques by flow scheduling policy assuming λ = 1, µ = 50 and Cogent [13] topology over 500 time units.

  1 min ; P N +1 ← The minimum weight path from s N +1 to t N +1 with at most K -1 hops (if multiple minimum weight paths exist, choose the one with minimum hops);

NP-Hardness proof dropped for brevity. It can be shown that the Set Cover problem can be reduced to Problem

In particular, we are looking for a subset of flows γ ⊆ F F F with minimum total remaining data units where there exists a path from s N +1 to t N +1 removing all edges that have a flow in F F F -γ.