EN FR
EN FR


Section: Scientific Foundations

Fast parametric estimation and its applications

Parametric estimation may often be formalized as follows:

y=F(x,Θ)+n,(1)

where:

  • the measured signal y is a functional F of the "true" signal x, which depends on a set Θ of parameters,

  • n is a noise corrupting the observation.

Finding a "good" approximation of the components of Θ has been the subject of a huge literature in various fields of applied mathematics. Most of those researches have been done in a probabilistic setting, which necessitates a good knowledge of the statistical properties of n. Our project is devoted to a new standpoint, which does not require this knowledge and which is based on the following tools, which are of algebraic flavor:

  • differential algebra (Differential algebra was introduced in nonlinear control theory by one of us almost twenty years ago for understanding some specific questions like input-output inversion. It allowed us to recast the whole of nonlinear control into a more realistic light. The best example is of course the discovery of flat systems, which are now quite popular in industry.), which plays with respect to differential equations a similar role that the commutative algebra plays with respect to algebraic equations;

  • module theory, i.e. linear algebra over rings, which are not necessarily commutative;

  • operational calculus, which is the most classical tool among control and mechanical engineers (Operational calculus is often formalized via the Laplace transform whereas the Fourier transform is today the cornerstone in estimation. Note that the one-sided Laplace transform is causal, but the Fourier transform over R is not.).

Linear identifiability

In the most problems, which appear in linear control as well as in signal processing, the unknown parameters are linearly identifiable: standard elimination procedures are yielding the following matrix equation

Pθ 1 θ r =Q,(2)

where:

  • θ i , 1ir represents unknown parameter,

  • P is a r×r square matrix and Q is a r×1 column matrix,

  • the entries of P and Q are finite linear combinations of terms of the form t ν d μ ξ dt μ , μ,ν0, where ξ is an input or output signal,

  • the matrix P is generically invertible, i.e., det(P)0.

How to deal with perturbations and noises?

With noisy measurements equation (2 ) becomes:

Pθ 1 θ r =Q+R,(3)

where R is a r×1 column matrix, whose entries are finite linear combinations of terms of the form t ν d μ η dt μ ,μ,ν0, where η is a perturbation or a noise.

Structured perturbations

A perturbation π is said to be structured if, and only if, it can be annihilated by a linear differential operator of the form finite a k (t)d k dt k , where a k (t) is a rational function of t, i.e. finite a k (t)d k dt k π=0. Note that many classical perturbations, like a constant bias, are annihilated by such an operator. An unstructured noise cannot be annihilated by a non-zero differential operator.

By well-known properties of the non-commutative ring of differential operators, we can multiply both sides of equation (3 ) by a suitable differential operator Δ such that equation (3 ) becomes:

ΔPθ 1 θ r =ΔQ+R ' ,(4)

where the entries of the r×1 column matrix R ' are unstructured noises.

Attenuating unstructured noises

Unstructured noises are usually dealt with stochastic processes like white Gaussian noises. They are considered here as highly fluctuating phenomena, which may therefore be attenuated via low pass filters. Note that no precise knowledge of the statistical properties of the noises is required.

Comments

Although the previous noise attenuation procedure (It is reminiscent to that the most practitioners in electronics are doing.) may be fully explained via formula (4 ), its theoretical comparison (Let us stress again that many computer simulations and several laboratory experiments have been already successfully achieved and can be quite favorably compared with the existing techniques.) with today's literature (Especially in signal processing.) has yet to be done. It will require a complete resetting of the notions of noises and perturbations. Besides some connections with physics, it might lead to quite new "epistemological" issues [101] .

Some hints on the calculations

The time derivatives of the input and output signals appearing in equations (2 ), (3 ), (4 ) can be suppressed in the two following ways which might be combined:

  • integrate both sides of the equation a sufficient number of times,

  • take the convolution product of both sides by a suitable low pass filter.

The numerical values of the unknown parameters Θ=(θ 1 ,,θ r ) can be obtained by integrating both sides of the modified equation (4 ) during a very short time interval.

A first, very simple example

Let us illustrate on a very basic example, the grounding ideas of the algebraic approach. For this purpose consider the first order linear system:

y ˙(t)=ay(t)+u(t)+γ 0 ,(5)

where a is an unknown parameter to be identified and γ 0 is an unknown constant perturbation. With the notations of operational calculus and y 0 =y(0), equation (5 ) reads:

sy ^(s)=ay ^(s)+u ^(s)+y 0 +γ 0 s(6)

where y ^(s) represents the Laplace transform of y(t).

In order to eliminate the term γ 0 , multiply first the both hand-sides of this equation by s and next take their derivatives with respect to s:

d dsssy ^(s)=ay ^(s)+u ^(s)+y 0 +γ 0 s(7)

2sy ^(s)+s 2 y ^ ' (s)=asy ^ ' (s)+y ^(s)+su ^ ' (s)+u ^(s)+y 0 .(8)

Recall that y ^ ' (s)dy ^(s) ds corresponds to -ty(t). Assume y 0 =0 for simplicity of presentation (If y 0 0 one has to take above derivatives of order 2 with respect to s, in order to eliminate the initial condition.). Then for any ν>0,

s -ν 2sy ^(s)+s 2 y ^ ' (s)=s -ν a(sy ^ ' (s)+y ^(s))+su ^ ' (s)+u ^(s).(9)

For ν=3, we obtained the estimated value a:

a=2 0 T dλ 0 λ y(t)dt- 0 T ty(t)dt+ 0 T dλ 0 λ tu(t)dt- 0 T dλ 0 λ dσ 0 σ u(t)dt 0 T dλ 0 λ dσ 0 σ y(t)dt- 0 T dλ 0 λ ty(t)dt(10)

Since T>0 can be very small, estimation via (10 ) is very fast.

Note that equation (10 ) represents an on-line algorithm, which involves only two kinds of operations on u and y: (1) multiplications by t, and (2) integrations over a pre-selected time interval.

If we now consider an additional noise of zero mean in (5 ), say:

y ˙(t)=ay(t)+u(t)+γ 0 +n(t),(11)

it can be considered as a fast fluctuating signal. The order ν in (9 ) determines the order of iterations in the integrals (3 integrals in (10 )). Those iterated integrals are low-pass filters which are attenuating the fluctuations.

This example, even simple, clearly demonstrates how algebraic techniques proceed:

  • they are algebraic: operations on s-functions;

  • they are non-asymptotic: parameter a is obtained from (10 ) in a finite time;

  • they are deterministic: no knowledge of the statistical properties of the noise n is required.

A second simple example, with delay

Consider the first order, linear system with constant input delay (This example is taken from  [93] . For further details, we suggest the reader to refer to it.):

y ˙(t)+ay(t)=y(0)δ+γ 0 H+bu(t-τ).(12)

Here we use a distributional-like notation, where δ denotes the Dirac impulse and H is its integral, i.e. the Heaviside function (unit step) (In this document, for the sake of simplicity, we make an abuse of the language since we merge in a single notation the Heaviside function H and the integration operator. To be rigorous, the iterated integration (k times) corresponds, in the operational domain, to a division by s k , whereas the convolution with H (k times) corresponds to a division by s k /(k-1)!. For k=0, there is no difference and H*y realizes the integration of y. More generally, since we will always apply these operations to complete equations (left-and right-hand sides), the factor (k-1)! makes no difference.). Still for simplicity, we suppose that the parameter a is known. The parameter to be identified is now the delay τ. As previously, γ 0 is a constant perturbation, a, b, and τ are constant parameters. Consider also a step input u=u 0 H. A first order derivation yields:

y ¨+ay ˙=ϕ 0 +γ 0 δ+bu 0 δ τ ,(13)

where δ τ denotes the delayed Dirac impulse and ϕ 0 =(y ˙(0)+ay(0))δ+y(0)δ (1) , of order 1 and support {0}, contains the contributions of the initial conditions. According to Schwartz theorem, multiplication by a function α such that α(0)=α ' (0)=0, α(τ)=0 yields interesting simplifications. For instance, choosing α(t)=t 3 -τt 2 leads to the following equalities (to be understood in the distributional framework):

t 3 [y ¨+ay ˙]=τt 2 [y ¨+ay ˙],bu 0 t 3 δ τ =bu 0 τt 2 δ τ .(14)

The delay τ becomes available from k1 successive integrations (represented by the operator H), as follows:

τ=H k (w 0 +aw 3 ) H k (w 1 +aw 2 ),t>τ,(15)

where the w i are defined using the notation z i =t i y by:

w 0 =t 3 y (2) =-6z 1 +6z 2 (1) -z 3 (2) ,w 1 =t 2 y (2) =-2z 0 +4z 1 (1) -z 2 (2) ,w 2 =t 2 y (1) =2z 1 -z 2 (1) ,w 3 =t 3 y (1) =3z 2 -z 3 (1) .

These coefficients show that k2 integrations avoid any derivation in the delay identification.

Figure 1. Delay τ identification from algorithm (15 )
IMG/tau2.png

Figure 1 gives a numerical simulation with k=2 integrations and a=2,b=1,τ=0.6, y(0)=0.3,γ 0 =2,u 0 =1. Due to the non identifiability over (0,τ), the delay τ is set to zero until the numerator or the denominator in the right hand side of (15 ) reaches a significant nonzero value.

Again, note the realization algorithm (15 ) involves two kinds of operators: (1) integrations and (2) multiplications by t. It relies on the measurement of y and on the knowledge of a. If a is also unknown, the same approach can be utilized for a simultaneous identification of a and τ. The following relation is derived from (14 ):

τ(H k w 1 )+aτ(H k w 2 )-a(H k w 3 )=H k w 0 ,(16)

and a linear system with unknown parameters (τ,aτ,a) is obtained by using different integration orders:

H 2 w 1 H 2 w 2 H 2 w 3 H 3 w 1 H 3 w 2 H 3 w 3 H 4 w 1 H 4 w 2 H 4 w 3 τ ^aτ ^-a ^=H 2 w 0 H 3 w 0 H 4 w 0 .

The resulting numerical simulations are shown in Figure 2 . For identifiability reasons, the obtained linear system may be not consistent for t<τ.

Figure 2. Simultaneous identification of a and τ from algorithm (16 )
IMG/exemple1AetTau.png