# Interval analysis for guaranteed non-linear parameter and state estimation

Interval analysis for guaranteed non-linear parameter and state estimation Mathematical and Computer Modelling of Dynamical Systems Vol. 11, No. 2, June 2005, 171 – 181 Interval analysis for guaranteed non-linear parameter and state estimation MICHEL KIEFFER* AND ERIC WALTER Laboratoire des Signaux et Syste` mes – CNRS – Supe´ lec – Universite´ Paris-Sud Plateau de Moulon, F-91192 Gif-sur-Yvette Cedex, France This paper presents some tools based on interval analysis for guaranteed non-linear parameter and state estimation in a bounded-error context. These tools make it possible to compute outer (and sometimes inner) approximations of the set of all parameter or state vectors that are consistent with the model structure, measurements and noise bounds. Keywords: Bounded errors, identiﬁcation, interval analysis, non-linear estimation, parameter bounding, state bounding. 1. Introduction Parameter and state estimation problems are encountered when modelling processes that involve uncertain quantities to be estimated from measurements. In this paper, it is assumed that all uncertain quantities (state perturbation, measurement noise, parameter vector, etc.) are bounded and belong to known sets. Many tools are available for characterizing the set of all values of the parameter or state vector that are consistent with this hypothesis. Most of these methods apply when the model output is linear in the parameters or initial state. The set to be characterized is then usually enclosed in simple geometric shapes such as polytopes, parallelotopes or ellipsoids (see, e.g.  and references therein). When the model under study is non- linear, the situation is much more complicated, as the set to be characterized may be non-convex and may even consist of several disconnected components. As we shall see, interval analysis provides tools for enclosing such sets in unions of non- overlapping interval vectors, the precision of the approximation being tuned by the *Corresponding author. Email: kieﬀer@lss.supelec.fr; walter@lss.supelec.fr Such a situation may even occur in the linear case, when the perturbation or measurement-noise sets are not convex or consist of disconnected parts. Mathematical and Computer Modelling of Dynamical Systems ISSN 1387-3954 print/ISSN 1744-5051 online ª 2005 Taylor & Francis Group Ltd http://www.tandf.co.uk/journals DOI: 10.1080/13873950500068807 172 M. Kieﬀer and E. Walter user. The aim of this paper is to summarize some recent results on guaranteed parameter and state estimation in this non-linear bounded-error context. Interval analysis was initially developed to assess the numerical errors resulting from the use of ﬁnite-precision arithmetic . It now makes it possible to obtain guaranteed solutions to standard engineering problems such as computing all solutions of a non- linear system of equations , or all global minimizers of a cost function [4,5] or computing sets guaranteed to contain all solutions of sets of non-linear inequalities  or enclosing the solution of a system of ordinary diﬀerential equations [7 – 10]. Libraries for interval analysis developed in object-oriented languages [11,12] or in environments such as Matlab  have helped disseminate these techniques. Section 2 will explain how parameter and state estimation in a bounded-error context can be formulated as set-estimation problems. In Section 3, interval analysis will be presented and tools for obtaining outer approximations of the image and reciprocal image of a set by a given function will be introduced. Section 4 will be devoted to the application of these tools to the problems mentioned in Section 2. Simple illustrative examples and pointers to real-world applications will also be provided before drawing some conclusions and perspectives. 2. Parameter and state estimation as set estimation problems The aim of this section is to show that non-linear parameter and state estimation can be formulated in a bounded-error context as problems of characterization of sets. Obtaining approximate outer (and sometime inner) approximations of these sets is possible with the tools provided by interval analysis to be presented in Section 3. 2.1 Non-linear Parameter Estimation Consider a system with known input u(t) and output y(t). For the sake of simplicity of notation, the dependence of y on u will be omitted. Deﬁne the output error vðp; tÞ¼ yðtÞ y ðp; tÞ; ð1Þ where p is a vector of unknown constant parameters and the model output y (p, t) may be computed by a function or ﬁnite algorithm or as the solution of a diﬀerential equation. Assume that v(t) and vðtÞ are known lower and upper bounds for the acceptable output errors. Such bounds may, for instance, correspond to bounded measurement noise. Assume further that the system output has been measured at discrete time instants t ,. . .,t . Bounded-error parameter estimation then corresponds to characteriz- 1 N ing the set P ¼fpjvðt Þ yðt Þ y ðp; t Þ vðt Þ; i ¼ 1; ... ; Ng; ð2Þ i i i i which contains all values of the parameter vector p that are consistent with the measurements and acceptable error bounds, and where all inequalities have to be understood componentwise. P may be rewritten as N N \ \ P ¼fg pjy ðp; t Þ2½yðt Þ ¼ P ; ð3Þ m i i i i¼1 i¼1 Interval analysis for guaranteed non-linear estimation 173 where [y(t )] is the interval vector with bounds yðt Þ vðt Þ and y(t )7v(t ). Bounded- i i i i i error parameter estimation thus amounts to the characterization of the intersection of the sets P associated with each measurement as illustrated by Figure 1(a). 2.2 Robust Non-linear Parameter Estimation In some situations, illustrated in Figure 1(b), the intersection of all P s is empty, which means that the hypotheses are contradictory. This may be due, for example, to overoptimistic noise bounds or to sensor failures at given time instants. The corresponding measurements will be called outliers. To obtain a solution set despite the presence of outliers, one may switch to robust parameter estimation. Consider the following test associated with the set P in (3) 1if y ðp; tÞ2½yðt Þ i i tðp; tÞ¼ 0 else: P can be rewritten as () P ¼ p j tðpÞ¼ tðp; tÞ¼ N i¼1 To protect the estimator against at most n outliers, it suﬃces to characterize the set P ¼fp j tðpÞ N  ngð4Þ Therefore, robust parameter estimation in a bounded-error context can again be cast into the framework of set characterization. For the case illustrated by Figure 1(b), an estimate P robust to at most one outlier obtained according to (4) is represented in Figure 1(c). 2.3 Non-linear State Estimation When the parameter vector varies with time and an equation describing the possible variations is available, one may switch to non-linear state estimation. Here, for the sake Figure 1. Bounded-error parameter estimation: (a) solution set P; (b) when outliers are present, P may be empty; (c) new solution set P robust to one outlier 1 174 M. Kieﬀer and E. Walter of simplicity, the state-space model is supposed to consist of a discrete-time diﬀerence equation x ¼ f ðx ; w ; uÞð5Þ kþ1 k k k k with initial condition x and an observation equation y ¼ g ðx Þþ v ð6Þ k k k k where x is the state vector to be estimated at time t , f and g are known functions or k k k k ﬁnite algorithms, u is some known input and w and v represents some state k k k perturbations and measurement noise. State perturbations account for the fact that (5) is only an approximation of reality. Measurement noise describes the uncertainty in the output measurements. Assume that w and v are bounded with known bounds, i.e. w 2 [w ]and v 2 [v ]. k k k k k k The initial state vector x is also assumed to belong to some known set X and at each 0 0 discrete instant of time t , the output vector y becomes available. With these k k assumptions, causal bounded-error state estimation aims at ﬁnding the set X of all state vectors that are consistent with the model structure (5) and (6), the measurements and noise bounds up to time t . An idealized recursive algorithm can easily be derived to achieve this task For ‘=0 to k7 1 (a) Prediction step : compute X ¼ f ðX ;½w ; u Þ , ‘ ‘ ‘ ‘ (b) Correction step : compute X ={x 2X j g (x)2 [y ]}, ‘+1 ‘ ‘ ‘ where ½y ¼½y  v ; y  v  . ‘ ‘ ‘ ‘ The correction step can be viewed as a problem of parameter estimation. The prediction step requires the computation of the direct image of the set X 6 [w ] by the ‘ ‘ function f . 3. Interval analysis for guaranteed set computation In Section 2, it has been shown that non-linear parameter estimation, robust parameter estimation and state estimation may be formulated as set characterization problems. These sets may not have simple shapes and may consist of several disconnected components. In this section, after presenting basic tools of interval analysis, two algorithms will be introduced that can compute outer approximations of sets of arbitrary shape. 3.1 Basic Tools An interval ½x¼½x; x is a closed and connected subset of R; it may be characterized by its lower and upper bounds x and x or equivalently by its centre cð½xÞ ¼ ðx þ xÞ=2 and width wð½xÞ ¼ x  x: Arithmetical operations on intervals can be deﬁned by Interval analysis for guaranteed non-linear estimation 175 8 2 fþ;;;=g;½x ½y¼fx yjx 2½x; y 2½yg Obtaining an interval corresponding to ½x ½y is easy for the ﬁrst three operators as ½xþ½y¼½x þ y; x þ y; ½x½y¼ ½x;y; x  y; ½x½y¼½minðxy; xy; xy; xyÞ; maxðxy; x; y; xy; xyÞ For division, when 02= [y], ½x=½y¼ ½minðx=y; x=y; x=y; x=yÞ; maxðx=y; x=y; x=y; x=yÞ and extended intervals have to be introduced when 02 [y] see, e.g. . More generally, the interval counterpart of a real-valued function is an interval- valued function deﬁned as fð½xÞ ¼ ½ffðxÞj 2 ½xg where [S] is the interval hull of S, i.e. the smallest interval that contains it. Interval counterparts to continuous elementary functions are easily obtained. For monotonic functions only computations on bounds are required expð½xÞ ¼ ½expðxÞ; expðxÞ logð½xÞ ¼ ½logðxÞ; logðxÞ if x > 0 For non-monotonic elementary functions, such as the trigonometric functions, algorithmic deﬁnitions are still easily obtained. For instance, the interval square function can be deﬁned by 2 2 ½0; maxðx ; x Þ; if 0 2½x ½x ¼ 2 2 2 2 ½minðx ; x Þ; maxðx ; x Þ else For more complicated functions, it is usually no longer possible to evaluate their interval counterpart, hence the importance of the concept of inclusion function.An inclusion function [f](.) for a function f(.) deﬁned over a domain D R is such that the image of an interval by this function is an interval, guaranteed to contain the image of the same interval by the original function: 8½x D; fð½xÞ ½fð½xÞ: ð7Þ This inclusion function is convergent if lim w([f]([x])) = 0 and inclusion w([x])?0 monotonic if [x] [y] implies [f]([x]) [f]([y]). Various techniques are available for building convergent and inclusion-monotonic inclusion functions. Among them, the simplest is to replace all occurrences of the real variable by its interval counterpart which results in what is called a natural inclusion function. Example 1 Consider the function fðxÞ¼ x  3ðx  expðxÞÞ An inclusion function for f is 176 M. Kieﬀer and E. Walter ½f ð½xÞ ¼ ½x  3ð½x expð½xÞÞ Evaluate [f] over [0, 1], ½f ð½0; 1Þ ¼ ½0; 1  3ð½0; 1 expð½0; 1ÞÞ ¼½0; 1 3ð½0; 1½1; eÞ ¼ ½0; 1 3½e; 0 ¼½0; 1þ½0; 3e¼½0; 3e þ 1 ½0; 9:16 Compare with fð½0; 1Þ ¼ ½3;2 þ 3e ½3; 6:16 Of course, f([0,1]) [f]([0,1]). When the inclusion in (7) becomes an equality, the inclusion function is minimal. Usually, some pessimism is introduced by the inclusion function, as in Example 1. This pessimism is due to the fact that each occurrence of the interval variable is considered as independent from the others. Various approaches may be considered to reduce pessimism. A ﬁrst one is to reduce the number of occurrences of the variable by symbolic manipulations, as illustrated by the following example. Example 2 The function fðxÞ¼ x  3ðx  expðxÞÞ can be rewritten as fðxÞ¼ðx  3=2Þ  9=4 þ 3 expðxÞ A new natural inclusion function can then be derived as ½f ð½xÞ ¼ ð½x 3=2Þ  9=4 þ 3 expð½xÞ Then ½f ð½0; 1Þ ¼ ð½0; 1 3=2Þ  9=4 þ 3 expð½0; 1Þ ¼½1; 3e ½1; 8:16 and f([0,1]) [f] ([0,1]) [f] ([0,1]),so [f] is less pessimistic than [f] . 2 1 2 1 A second approach is to use more sophisticated inclusion functions, based, e.g. on Taylor expansions (see , ,  or ). A third approach resorts to the bisection of [x] into two smaller intervals [x ]and [x ] 1 2 and uses the property that for inclusion-monotonic inclusion functions fð½xÞ ½½fð½x Þ [ ½fð½x Þ ½fð½xÞ ð8Þ 1 2 This technique will be intensively used in the Sivia and ImageSp algorithms to be presented. Despite the fact that inclusion functions usually only provide an outer approxima- tion of the range of a function over an interval, this approximation may allow Interval analysis for guaranteed non-linear estimation 177 mathematical results to be proven numerically. It is for example possible to prove that a function f(x) does not vanish over a given interval [x] using one of its inclusion functions. If 02= [f]([x]), then as f([x]) [f]([x]), there will be no x2 [x] such that f(x)=0. Such a proof may even be established on a computer with ﬁnite-precision ﬂoating-point representation of real numbers, provided that every computation with intervals is outwards rounded, as illustrated by the following example. Example 3 Assume that all computations of Example 2 are performed on a computer that has a two-digit representation of numbers. Then the interval provided by an inclusion function evaluated with outwards rounding would be ½f ð½0; 1Þ ¼ ð½0; 1 3=2Þ  9=4 þ 3 expð½0; 1Þ ¼½1:5;0:5 ½2:2; 2:3þ 3½1:0; 2:8 ¼½0:2; 2:3½2:2; 2:3þ½3:0; 8:4¼ ½0:9; 8:5 The relation f([0,1]) [f] ([0,1]) is still true, and it remains possible to prove that f(x) does not vanish on [0,1]. More details on guaranteed computation with intervals using a ﬁnite-precision representation of numbers may be found in  or . An interval vector (or box)[x] may equivalently be seen as the Cartesian product of scalar intervals ½x¼ ½x  ½x 1 n or as a vector with interval components ½x¼ ð½x ; ... ;½x Þ 1 n The width w([x]) of a box is the maximum of the widths of its components. Inclusion functions for vector functions of interval vectors are easily deﬁned. Unions of non- overlapping boxes will be called subpavings. 3.2 Inverse Image Evaluation (Set Inversion) The Sivia algorithm (for Set Inverter Via Interval Analysis ) is dedicated to the characterization of sets deﬁned by P ¼fpjfðpÞ 0gð9Þ where f(.) is any function for which an inclusion function [f](.) is available. An initial search box [p ] to which the search will be restricted must also be provided. Sivia is indeed a set inverter as it characterizes P ¼ f ð½0;þ1½ Þ It recursively builds two unions of non-overlapping boxes (subpavings) P and P such that P 178 M. Kieﬀer and E. Walter according to the following algorithm: 1. Initialization: put [p ] in the list L of boxes to be treated. 2. Take a box [p] out of L. 6 n a. If [f]([p]) [0, +?[ , store [p]in P and in P. 6 n b. If [f]([p]) \ ([0, +?[ =;,goto3. c. If w([p])5 e, store [p]in P else bisect [p] into [p ] and [p ] and store them in L. 1 2 3. If L = ; go to 2. To ensure termination after a ﬁnite number of iterations, a box [p] may be bisected only if its width is larger than some prespeciﬁed precision parameter e. 3.3 Direct Image Evaluation Direct image evaluation is the characterization of Y ¼ fðXÞ when X is a known set and f a known function. This task is usually quite complex. The algorithm ImageSp, again based on interval analysis, builds a subpaving Y such that Y when X is itself a subpaving and an inclusion function [f] is available for f. ImageSp consists of three steps: 1. Mincing: all boxes in X are bisected in order to obtain a subpaving X’ consisting only of boxes of width less than a prespeciﬁed precision parameter e. 2. Evaluation: the images of all boxes of X’ are evaluated using [f] and stored into a list of image boxes Y. 3. Regularization: all boxes in Y are merged into a new subpaving Y. Y is guaranteed to contain Y. The precision of the approximation is controlled by e. For more details, see  or . 4. Applications Three very simple illustrative examples will be presented here. Pointers to real-life estimation problems solved using interval computation will be also provided. 4.1 Non-linear Parameter Estimation The deﬁnition of the set to be characterized in a parameter estimation problem (2) is easily transformed into an expression similar to (9). Sivia is thus a natural tool for guaranteed non-linear parameter estimation in a bounded-error context. Consider a system on which the data of Table 1 have been collected. This system is modelled by a sum of two exponentials y ðp; tÞ¼ p expðp tÞþ p expðp tÞ; m 1 2 3 4 Interval analysis for guaranteed non-linear estimation 179 and the output error is considered acceptable if it belongs to [7 0.01, 0.01] (corresponding to the accuracy of the sensor). 6 4 The Sivia algorithm is applied with [p ] = [0, 5] and e = 0.005. The natural inclusion function is employed for y (p, t). The projections onto the (p , p ) and (p , m 1 2 3 p ) planes of an outer approximation of the solution obtained in less than 7 s on an Athlon at 1 GHz are presented on Figure 2. This solution consists of two disconnected components, as expected, since (p , p )and (p , p ) can be exchanged without 1 2 3 4 modifying the behaviour of the model. Applications to real-life problems may be found, e.g. in  for the estimation of electrochemical parameters and in  for the estimation of thermal characteristics. 4.2 Robust Non-linear Estimation Assume now that a sensor failure forces y(2) to 0. If the Sivia algorithm is applied again with this new measurement set, the solution is proved to be empty in less than 0.1 s. For robust parameter estimation, the deﬁnition of the set (4) to be characterized can again be expressed as in (9). Sivia can thus again be employed for the characterization of P . Using the same tuning of the algorithm as before, with n = 1, the projection on the (p , p ) and (p , p ) planes of an outer approximation of the solution obtained in 1 2 3 4 less than 8 s is represented in Figure 3. The volume of the set obtained is slightly larger Table 1. Data measured on the system of Section 4.1 t(s) 01234 5 y (t) 3.01 1.56 0.93 0.65 0.49 0.39 Figure 2. Projections on the (p , p ) and (p , p ) planes of an outer approximation of the solution set for the 1 2 3 4 parameter estimation problem (no outlier) Figure 3. Projections on the (p , p ) and (p , p ) planes of an outer approximation of the solution set for the 1 2 3 4 robust parameter estimation problem (one outlier) 180 M. Kieﬀer and E. Walter than in Section 4.1. This is due to the fact that less information has been employed for the characterization of P , as one datum had to be discarded to get a non-empty solution set. Applications of robust parameter estimation based on interval analysis may be found, e.g. in  in the context of robot localization in the presence of unreliable ultrasonic measurements. 4.3 Non-linear State Estimation The implementable version of the recursive non-linear state estimation algorithm presented in Section 2.3 consists of two steps. The prediction step involves the computation of an outer approximation of the direct image of a set, which can be achieved by the ImageSp algorithm. The correction step can be viewed as a parameter estimation problem, which is solved using Sivia. For more details, see . Application to real-life problems may be found, e.g. in  for robot tracking and in  in conjunction with guaranteed ODE solvers for the estimation of the state of a waste-water treatment unit. 5. Conclusions Guaranteed techniques for non-linear parameter and state estimation in a bounded- error context have been presented. The main tools for this purpose are guaranteed algorithms for inverse and direct image evaluation. These algorithms are based on interval analysis, which we believe to be an extremely promising approach for the investigation of the properties of non-linear models in a bounded-error context. Global optimization algorithms based on interval analysis [4,14] can also be used to ﬁnd optimal estimates in a context where uncertainty is characterized probabilistically, for instance by allowing a guaranteed computation of maximum-likelihood estimates. Pointers to software and many details about implementation issues may be found elsewhere, e.g. in  or on the web at http://www.cs.utep.edu/interval-comp/ main.html. References  Milanese, M., Norton, J., Piet-Lahanier, H. and Walter, E. (Eds), 1996, Bounding Approaches to System Identiﬁcation (New York: Plenum Press).  Moore, R.E., 1959, Automatic error analysis in digital computation. Technical Report LMSD-48421, Lockheed Missiles and Space Co., Palo Alto, CA.  Neumaier, A., 1990, Interval Methods for Systems of Equations (Cambridge: Cambridge University Press).  Hansen, E.R., 1992, Global Optimization Using Interval Analysis (New York: Marcel Dekker).  Ratschek, H. and Rokne, J., 1988, New Computer Methods for Global Optimization (Chichester: Ellis Horwood).  Jaulin, L. and Walter, E., 1993, Automatica, 29, 1053 – 1064.  Berz, M. and Makino, K., 1998, Reliable Computing, 4, 361 – 369.  Lohner, R., 1992, In J.R. Cash and I. Gladwell (Eds) Computational Ordinary Diﬀerential Equations, pages 425 – 435 (Oxford: Clarendon Press).  Moore, R.E., 1979, Methods and Applications of Interval Analysis (Philadelphia: SIAM).  Nedialkov, N.S. and Jackson, K.R., 2001, In U. Kulisch, R. Lohner, and A. Facius (Eds) Perspectives on Enclosure Methods, pages 219 – 264 (Vienna: Springer-Verlag).  Kearfott, R.B., 1995, ACM Transactions on Mathematical Software, 21, 63 – 78. Interval analysis for guaranteed non-linear estimation 181  Klatte, R., Kulisch, U., Wieth, A., Lawo, C. and Rauch, M., 1993, C-xsc: A C++ A Class Library for Extended Scientiﬁc Computing (Berlin: Springer-Verlag).  Rump, S.M., 2001, In J.*Grabmeier, E. Kaltofen, and V. Weispfennig (Eds) Handbook of Computer Algebra: Foundations, Applications, Systems (Heidelberg: Springer-Verlag).  Jaulin, L., Kieﬀer, M., Didrit, O. and Walter, E., 2001, Applied Interval Analysis (London: Springer- Verlag).  Ratschek, H. and Rokne, J., 1984, Computer Methods for the Range of Functions (Chichester: Ellis Horwood).  Chiriaev, D. and Walster, G.W., 1998, Interval arithmetic speciﬁcation. Available online at: http:// www.mscs.edu/globsol/walster-papers.html.  Lerch, M. and Wolﬀ von Gudenberg, J., 2000, ﬁ_lib++: Speciﬁcation, implementation and test of a library for extended interval arithmetic. In Proc. 4th Conference on Real Numbers and Computers, Dagstuhl, Germany. Available online at: http://www.imada.sdu.dk/\symbol{126}kornerup/RNC4/ papers/p03.ps.  Kieﬀer, M., Jaulin, L. and Walter, E., 2002, International Journal of Adaptative Control and Signal Processing, 6, 193 – 218.  Braems, I., Berthier, F., Jaulin, L., Kieﬀer, M. and Walter, E., 2001, Journal of Electroanalytical Chemistry, 495,1–9.  Braems, I., Jaulin, L., Kieﬀer, M., Ramdani, N. and Kieﬀer, M., 2003, Reliable parameter estimation in presence of uncertain variables that are not estimated. In Proceedings of the 13th IFAC Symposium on System Identiﬁcation SYSID 2003, Rotterdam, The Netherlands, Elsevier Science Ltd, Oxford, pages 1856 – 1861.  Kieﬀer, M., Jaulin, L., Walter, E. and Meizel, D., 2000, Reliable Computing, 6, 337 – 362.  Kieﬀer, M., Jaulin, L., Walter, E. and Meizel, D., 1999, Guaranteed mobile robot tracking using interval analysis. In Proc. MISC’99 Workshop on Applications of Interval Analysis to Systems and Control, pages 347 – 359, Girona, Spain.  Kieﬀer, M., and Walter, E. 2004, Guaranteed nonlinear state estimator for cooperative systems. Numerical Algorithms, 37, 187 – 198. http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png Mathematical and Computer Modelling of Dynamical Systems Taylor & Francis

# Interval analysis for guaranteed non-linear parameter and state estimation

, Volume 11 (2): 11 – Jun 1, 2005
11 pages      