Particle observers for contracting dynamical systems

07/11/2017
Publication GSI2017
OAI : oai:www.see.asso.fr:17410:22540
contenu protégé  Document accessible sous conditions - vous devez vous connecter ou vous enregistrer pour accéder à ou acquérir ce document.
- Accès libre pour les ayants-droit
 

Résumé

In the present paper we consider a class of partially observed dynamical systems. As in the Rao-Blackwellized particle lter (RBPF) paradigm (see e.g., [Doucet et al. (2000)]), we assume the state x can be broken into two sets of variables x = (z; r) and has the property that conditionally on z the system's dynamics possess geometrical contraction properties, or is amenable to such a system by using a nonlinear observer whose dynamics possess contraction properties. Inspired by the RBPF we propose to use particles to approximate the r variable and to use a simple copy of the dynamics (or an observer) to estimate the rest of the state. 
This has the bene ts of 1- reducing the computational burden (a particle filter would sample the variable x also), which is akin to the interest of the RBPF, 2- coming with some indication of stability stemming from contraction (actual proofs of stability seem difficult), and 3- the obtained filter is well suited to systems where the dynamics of x conditionally on z is precisely known and the dynamics governing the evolution of z is quite uncertain.

Particle observers for contracting dynamical systems

Collection

application/pdf Particle observers for contracting dynamical systems Silvère Bonnabel, Jean-Jacques Slotine
Détails de l'article
contenu protégé  Document accessible sous conditions - vous devez vous connecter ou vous enregistrer pour accéder à ou acquérir ce document.
- Accès libre pour les ayants-droit

Particle observers for contracting dynamical systems

Média

Voir la vidéo

Métriques

0
0
334.2 Ko
 application/pdf
bitcache://1a8fe37610efc0479801cf7ba126f9dd25f0073f

Licence

Creative Commons Aucune (Tous droits réservés)

Sponsors

Sponsors Platine

alanturinginstitutelogo.png
logothales.jpg

Sponsors Bronze

logo_enac-bleuok.jpg
imag150x185_couleur_rvb.jpg

Sponsors scientifique

logo_smf_cmjn.gif

Sponsors

smai.png
logo_gdr-mia.png
gdr_geosto_logo.png
gdr-isis.png
logo-minesparistech.jpg
logo_x.jpeg
springer-logo.png
logo-psl.png

Organisateurs

logo_see.gif
<resource  xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance"
                xmlns="http://datacite.org/schema/kernel-4"
                xsi:schemaLocation="http://datacite.org/schema/kernel-4 http://schema.datacite.org/meta/kernel-4/metadata.xsd">
        <identifier identifierType="DOI">10.23723/17410/22540</identifier><creators><creator><creatorName>Silvère Bonnabel</creatorName></creator><creator><creatorName>Jean-Jacques Slotine</creatorName></creator></creators><titles>
            <title>Particle observers for contracting dynamical systems</title></titles>
        <publisher>SEE</publisher>
        <publicationYear>2018</publicationYear>
        <resourceType resourceTypeGeneral="Text">Text</resourceType><dates>
	    <date dateType="Created">Wed 7 Mar 2018</date>
	    <date dateType="Updated">Wed 7 Mar 2018</date>
            <date dateType="Submitted">Fri 17 Aug 2018</date>
	</dates>
        <alternateIdentifiers>
	    <alternateIdentifier alternateIdentifierType="bitstream">1a8fe37610efc0479801cf7ba126f9dd25f0073f</alternateIdentifier>
	</alternateIdentifiers>
        <formats>
	    <format>application/pdf</format>
	</formats>
	<version>37251</version>
        <descriptions>
            <description descriptionType="Abstract">In the present paper we consider a class of partially observed dynamical systems. As in the Rao-Blackwellized particle lter (RBPF) paradigm (see e.g., [Doucet et al. (2000)]), we assume the state x can be broken into two sets of variables x = (z; r) and has the property that conditionally on z the system's dynamics possess geometrical contraction properties, or is amenable to such a system by using a nonlinear observer whose dynamics possess contraction properties. Inspired by the RBPF we propose to use particles to approximate the r variable and to use a simple copy of the dynamics (or an observer) to estimate the rest of the state. <br />
This has the bene ts of 1- reducing the computational burden (a particle filter would sample the variable x also), which is akin to the interest of the RBPF, 2- coming with some indication of stability stemming from contraction (actual proofs of stability seem difficult), and 3- the obtained filter is well suited to systems where the dynamics of x conditionally on z is precisely known and the dynamics governing the evolution of z is quite uncertain.
</description>
        </descriptions>
    </resource>
.

Particle observers for contracting dynamical systems Silvère Bonnabel1 and Jean-Jacques Slotine2 1 MINES ParisTech, PSL Reasearch University, Centre for Robotics, 60 bd Saint-Michel, 75006 Paris, France, silvere.bonnabel@mines-paristech.fr, 2 Department of Mechanical Engineering, Massachusetts Institute of Technology, MA, 02139 USA, jjs@mit.edu. Abstract. In the present paper we consider a class of partially observed dynamical systems. As in the Rao-Blackwellized particle filter (RBPF) paradigm (see e.g., [Doucet et al. (2000)]), we assume the state x can be broken into two sets of variables x = (z, r) and has the property that conditionally on z the system’s dynamics possess geometrical contraction properties, or is amenable to such a system by using a nonlinear observer whose dynamics possess contraction properties. Inspired by the RBPF we propose to use particles to approximate the r variable and to use a simple copy of the dynamics (or an observer) to estimate the rest of the state. This has the benefits of 1- reducing the computational burden (a particle filter would sample the variable x also), which is akin to the interest of the RBPF, 2- coming with some indication of stability stemming from contraction (actual proofs of stability seem difficult), and 3- the obtained filter is well suited to systems where the dynamics of x conditionally on z is precisely known and the dynamics governing the evolution of z is quite uncertain. 1 A primer on contraction theory 1.1 Background on contraction theory Consider a Riemannian manifold (M, g), where g denotes the metric. Consider local coordinates. In the present paper, we will simplify the exposure by system- atically assuming that M = Rn . The squared infinitesimal length is given by the quadratic form: kdxk2 = X 1≤i,j≤n gij(x)dxidxj The matrix G = (gij)1≤i,j≤n is called the Riemannian metric tensor and it generally depends on x. Now, consider the continuous time deterministic system described by the following ordinary differential equation (ODE) on Rn : d dt x = f(x), (1) with f a smooth nonlinear function satisfying the usual conditions for global ex- istence and unicity of the solution. For a detailed proof of the following theorem, see e.g., [Pham and Slotine (2013)]. Theorem 1 ([Lohmiller and Slotine(1998)]) Let Jf (x) denote the Jacobian matrix of f(x). Assume that M(x) = GT (x)G(x) is uniformly positive definite, and that G(x)Jf (x)G−1 (x) is uniformly negative definite, then all trajectories exponentially converge to a single trajectory. Moreover, the convergence rate is equal to λ > 0 which is the supremum over x of the largest eigenvalue of G(x)Jf (x)G−1 (x). More precisely, if a(t) and b(t) are two trajectories of (1), we have: dg(a(t), b(t)) ≤ dg(a(0), b(0))e−2λt , where dg denotes the Riemannian distance associated to metric g. 1.2 Nonlinear observers for contracting systems Consider the system (1) where x(t) ∈ RN , with partial observations y(t) = h(x(t)), (2) The goal of observer design, is to estimate in real time the unknown quantity x(t) with the greatest possible accuracy given all the measurements up to current time t. Assume that for a class of functions y(t), the dynamics d dt z = f(z) + K(z, y)(y − h(z)) (3) can be proved to be contractive with rate λ > 0. Then, the observer for the system (1)-(2) defined by d dt x̂ = f(x̂) + K(x̂, y)(y − h(x̂)), (4) possesses convergence properties. Indeed, as the simulated x̂(t) and the true trajectory x(t) are both solutions of equation (3), Theorem 1 applies and we have: dg(x̂(t), x(t)) ≤ dg(x̂(0), x(0))e−2λt . 2 The basic particle observer 2.1 The Rao-Blackwellized particle filter (RBPF) Consider a (discrete) Markov process rt of initial distribution p(r0) and tran- sition equation p(rt | rt−1). The variable rt is hidden, and assume we have as observation a random variable yt at time t, which is correlated with rt. The observations are assumed to be conditionally independent given the process rt. The goal of discrete time filtering is to infer online the hidden variables from the observations, that is, to compute: p(rt | y1:t), where y1:t = {y1, · · · , yt}, or more generally p(r1:t | y1:t). Assume now, that we also want to infer another related process zt, such that p(zt | y1:t, r1:t) can be analytically evaluated. This is typically the case using a Kalman filter when conditionally on r the system is linear and Gaussian. A simple version of the RBPF is given by Algorithm 1. Algorithm 1 RBPF with prior sampling (see e.g., [Doucet et al. (2000)]) Draw N particles from the prior initial distribution p(r0) loop Sample from the prior r (i) t ∼ p(rt | r (i) t−1), and let r (i) 1:t =  r (i) t , r (i) 1:t−1  Evaluate and update weights w (i) t = p(yt | y1:t−1, r (i) 1:t) w (i) t−1 Normalize weights w̃ (i) t = w (i) t [ P j w (j) t ]−1 The estimate of the expected value E (F(z, r)) of any function F is X i w̃ (i) t Ep(zt|y1:t,r (i) 1:t)  F(r (i) t , zt)  Resample if necessary, i.e., duplicate and suppress particles to obtain N random samples with equal weights (i.e., equal to 1/N). end loop 2.2 The particle observer for conditonnally contracting systems Consider a noisy dynamical system of the form d dt z = f(z, r) (5) d dt r = g(r) + w(t) (6) where f, g are smooth maps, w(t) is a process noise, and we have an initial prior distribution π0(z, r) at time t = 0. Assume one has access to discrete time uncertain measurements yn = h(ztn , rtn ) + Vn at times t0 < t1 < t2 < · · · , and where Vn are unknown independent identically distributed random variables with known density l, that is, p(y | z, r) = l(y − h(z, r)). We introduce the following definition. Definition 1. The system (5)-(6) is said to be a contraction conditionally on z if equation (5) is a contraction when r(t) is considered as a known input. The rationale of our particle observer is as follows. If r(t) were known, then, all trajectories of the system (5) would converge to each other due to the conditional contraction properties we assume. Thus, if we call ẑ(t) a solution of (5) associated to some trajectory {r(t)}t≥0, then asymptotically we have p(z(t) | {r(s)}0≤s≤t) ≈ δ(z(t) − ẑ(t)), which means that contrarily to the RBPF paradigm we can not compute the conditional distributions in closed form but we have access to relevant approximations to them. Thus, letting (ẑ (i) t , r (i) t ) be a solution to the stochastic differential equations (5)-(6), we have the following approximations that stem from the partial contraction properties of the system: p(zt | y1:t, r (i) 1:tn ) ≈ δ(zt − ẑ (i) t ), p(ytn | y1:tn−1 , r (i) 1:tn ) ≈ l(ytn − h(ẑ (i) tn , r (i) tn )). Thus, resorting to those approximation, and applying the RBPF methodology to the above system (5)-(6) we propose the following Algorithm 2. Algorithm 2 The PO with prior sampling Draw N particles (z (1) t0 , r (1) t0 ), · · · , (z (N) t0 , r (N) t0 ) from the prior initial distribution π0(z, r) loop Sample (z (i) tn , r (i) tn ) from the prior by numerically integrating the stochastic differential equations (5)-(6) from time tn−1 to tn. Evaluate and update weights w (i) t = l(ytn − h(z (i) tn , r (i) tn )) w (i) t−1. Numerically enforce that at least one weight is not equal to zero (i.e., if all weights are zero, set e.g. w (1) t = 1). Normalize weights w̃ (i) t = w (i) t [ P j w (j) t ]−1 . The estimate of the expected value E (F(z, r)) of any function F is approximated by X i w̃ (i) t F(r (i) t , z (i) t ) In particular the state is approximated by X i w̃ (i) t (r (i) t , z (i) t ) Resample if necessary, i.e., duplicate and suppress particles to obtain N random samples with equal weights (i.e., equal to 1/N). end loop 2.3 Some comments on the choice of model The relevance of system (5)-(6) is debatable, for the two following reasons. First, it might be surprising that the dynamics of z conditionally on r be deterministic, whereas the dynamics of r be noisy. Second, because it is rare to find systems that are naturally (conditionally) contracting. Both issues will be partly addressed in the extensions outlined in the sequel. At this stage, we can make the following comments regarding the first issue. Assume both equations (5)-(6) to be noisy. Then, thanks to the contraction property, the asymptotic distribution of z(t) conditionally on r(t) is not very dispersed if the process noise is moderate, see [Pham et al. (2009)]. So the method may yield good results in practice. Assume on the other hand, both equations (5)-(6) to be deterministic. Then, it is hopeless to estimate and track efficiently the state with a (RB) particle filter, as the state space will not be explored adequately. Indeed, because multiple copies are produced after each resampling step, the diversity of the particle system decreases to a few points, which can be very different from the true state. To solve this degeneracy problem, the regularized particle filter was proposed in [Musso and Oujdane (1998)]. Albeit debatable, this technique may yield good results in practice. Following this route, we can postulate noisy equation (6) to implement our particle filter. Remark 1 Note that, here we do not deal with parameter identification, as in e.g., [Saccomani et. al (2003)]. Although this might look similar, r(t) is not a parameter, preventing us to directly apply the results of e.g., [Wills et. al (2008)] 3 A chemical reactor example 3.1 Retained model Consider the exothermic chemical reactor of [Adebekun and Schork(1989)]. It was shown in [Lohmiller and Slotine(1998)] that, if the temperature T is known, and thus can be considered as an input, then the system is a contraction. But to achieve best performance, and filter the noise out of the temperature measure- ments, the temperature should be considered as a (measured) part of the state as in [Adebekun and Schork(1989)]. This leads to a system that is not a contrac- tion. To make our point, we even propose to slightly modify the temperature dynamics to make it clearly unstable, yielding the more challenging following system: d dt I = q(t) V (If − I) − kde− Ed RT (t) I (7) d dt M = q(t) V (Mf − M) − 2kpe− Ed RT (t) M2 I (8) d dt P = q(t) V (Pf − P) + kpe− Ed RT (t) M2 I (9) d dt T = βT + σ2w(t) (10) where w(t) is a white Gaussian standard noise, and σ2 > 0 a parameter encoding the noise amplitude. Letting Vn ∼ N(0, 1) a random standard centered Gaussian, we assume discrete temperature measurements of the form: yn = T(tn) + σ1Vn. (11) [Lohmiller and Slotine(1998)] already proved the system is contracting condi- tionally on T(t). Thus, we can use the method described in Algorithm 2. 3.2 Simulation results The true system is simulated according to the equations (7)-(8)-(9)-(10) where we turned the noise off in equation (10) (this means we started from a noise-free system for which the RBPF would not work properly, and used the regularization technique discussed at Section 2.3). The noisy output (11) was also simulated, where an observation is made every 5 steps. We chose σ̃2 = 0.1. Density l is dictated by the observation noise, that is, l(u) = 1 σ1 √ 2π exp(− u2 2σ2 1 ) with σ1 = 1K. To apply our methodology, we assume that we have plausible physical upper bounds on the concentrations inside, and denote them by Imax, Mmax, Pmax and we let π0 be the uniform distribution on the hyperrectangle [0, Imax]×[0, Mmax]× [0, Pmax] with a Dirac on the measured initial temperature. In the simulation, all those upper bounds are set equal to 4mol. q(t) V and T(t) are slowly oscillating around 1min−1 and 300K, we have kde+ Ed RT (t) ≈ 0.8min−1 and kpe+ Ed RT (t) ≈ 0.2L mol−1 min−1 . We also let β = 0.01min−1 . N=15 particles are used (which results in a very cheap to implement particle filter, as each particle is associated only to a naive observer). We resampled3 each time the number of effective particles 1/( PM 1 w2 j ) drops below N/4, i.e., 25% of the total population. The resampling step is necessary, so that all particles grad- ually improve their estimation of the temperature, allowing the concentrations to be well estimated in turn. The noise is efficiently filtered and all values asymptotically very well re- covered, although a very reduced number of particles is used (15 observers are running in parallel) and measured temperature is noisy. See Figures 1 and 2. 4 Conclusion We have proposed a novel method to estimate the state of a class of dynamical systems that possess partial contraction properties. The method has successfully been applied to a chemical reactor example. Possible extensions are twofold. First, if equation (5) is noisy, one can use the same RBPO. Using the result of [Pham et al. (2009),Pham and Slotine (2013)], we can have an approximation of the asymptotic variance associated to the distribution p(z(t) | {r(s)}0≤s≤t). Thus, a Gaussian approximation to this distribution can be leveraged to imple- ment a RBPF. Second, if f is not contracting conditionally on r, but, is amenable to it using an observer of the form d dt z = f(z) + K(z, y, r)(y − h(z)), then the method may still be applied. 3 i.e., draw N particles from the current particle set with probabilities proportional to their weights; replace the current particle set with this new one. Instead of setting the weights of the new particles equal to 1/N as in the standard methodology, we preferred in the simulations to assign them their former weight and then normalize. Fig. 1. Left: True concentrations (dashed lines) and trajectories of the 15 particles. We see the effect of resampling, that refocuses the bundle of trajectories on the fittest ones, when too many become unlikely. Right: True concentrations (dashed lines), and estimates of the particle observer (solid lines). The ideas introduced in this short paper might also be applied to differen- tially positive systems [Forni and Sepulchre(2016),Bonnabel et al.(2011)]. In the future, we would also like to study the behavior of particle filters for systems with contraction properties. A starting point could be to seek how to use the recent re- sults of [Pham et al. (2009),Pham and Slotine (2013),Tabareau et. al (2010)] on stochastic contraction. References [Adebekun and Schork(1989)] Adebekun, D.K. and Schork, F.J. (1989). Continuous solution polymerization reactor control. 1. nonlinear reference control of methyl methacrylate polymerization. Industrial & engineering chemistry research, 28(9), 1308–1324. [Aminzarey and Sontag(2014)] Aminzarey, Z. and Sontag, E.D. (2014). Contraction methods for nonlinear systems: A brief introduction and some open problems. In 53rd IEEE Conference on Decision and Control, 3835–3847. IEEE. [Bonnabel et al.(2011)] Bonnabel, S., Astolfi, A., and Sepulchre, R. (2011). Contrac- tion and observer design on cones. In IEEE Conference on Decision and Control 2011 (CDC11). [Doucet et al. (2000)] Doucet, A. , De Freitas, N. and Murphy, K. (2000). Rao- Blackwellised particle filtering for dynamic Bayesian networks. In Proceedings of the Sixteenth conference on Uncertainty in artificial intelligence, 176–183. [Forni and Sepulchre(2016)] Forni, F. and Sepulchre, R. (2016). Differentially positive systems. IEEE Transactions on Automatic Control, 61(2), 346–359. [Lohmiller and Slotine(1998)] Lohmiller, W. and Slotine, J. (1998). On contraction analysis for nonlinear systems. Automatica, 34(6), 683–696. [Musso and Oujdane (1998)] C. Musso and N. Oudjane (1998). Regularization schemes for branching particle systems as a numerical solving method of the nonlinear filtering problem. Proceedings of the Irish Signals and Systems Conference. Fig. 2. True (solid line), measured (noisy line), estimated (dashed line, output by the RBPO) temperatures. [Pham et al. (2009)] Pham, Q.C., Tabareau, N., and Slotine, J.J. (2009). A contrac- tion theory approach to stochastic incremental stability. IEEE Transactions on Automatic Control, 54(4), 816–820. [Pham and Slotine (2013)] Pham, Q. C. and Slotine, J. J. (2013). Stochastic contrac- tion in Riemannian metrics. arXiv preprint arXiv:1304.0340.. [Saccomani et. al (2003)] Saccomani, M. P., Audoly, S. and D’Angi, L. (2003). Param- eter identifiability of nonlinear systems: the role of initial conditions. Automatica, 39(4), 619-632. [Tabareau et. al (2010)] Tabareau, N., Slotine, J. J. and Pham, Q. C. (2010). How synchronization protects from noise. PLoS Comput Biol, 6(1), e1000637. [Wang and Slotine(2005)] Wang, W. and Slotine, J.J.E. (2005). On partial contraction analysis for coupled nonlinear oscillators. Biological cybernetics, 92(1), 38–53. [Wills et. al (2008)] Wills, A., Schn, T. B and Ninness, B. (2008). Parameter estima- tion for discrete-time nonlinear systems using EM. IFAC Proceedings Volumes, 41(2), 4012-4017.