As of the time of this writing, the goal of my work at Geostat has been, as a preliminary step, to reproduce the work of Almpanidis in voice activity detection or phone segementation using the DISTBIC-Γ method. For a review of this method, please see my previous works.

In this report I will describe the methods implemented or which were attempted in working towards a reproduction of the work of Almpanidis.

The DISTBIC-Γ method involves several steps. Most of these steps were implemented early in this project in C++, then translated into Matlab, which was a quick process. However, one step -the estimation of the parameters of the Generalized Γ distribution, has proven quite difficult to master and consequently has consumed most of my work. Because the other steps have already been described in my previous reports almost exactly as they are currently implemented, I will focus in this report on the methods used to estimate the GΓD.

Gradient Descent Algorithm

Theoretical Overview

This procedure is described in (Almpanidis, 2006(Almpanidis, , 2008) ) and (Shin, Chang, 2005). etc. as the method used by these researchers to estimate the parameters η β and γ for the GΓD.

Assume data are mutually independent Calculate the following statistics in order to estimate γ

S 1 (i) = (1 -ζ)S 1 (i -1) + ζ|x i | γ(i)
(1)

S 2 (i) = (1 -ζ)S 2 (i -1) + ζ log |x i | γ(i) (2) S 3 (i) = (1 -ζ)S 2 (i -1) + ζ|x i | γ(i) log |x i | γ(i) (3)
where the statistics in question are given by the following:

S 1 = (1/N ) N i=1 |x i | γ (4)
S 2 = (1/N ) N i=1 log |x i | γ (5)
S 3 = (1/N ) N i=1 |x i | γ log |x i | γ (6)
update γ as:

γ(i + 1) = γ(i) + µ 1 η(i) + S 2 (i) - S 3 (n) S 1 (i) (7)
Where ζ is a forgetting factor and µ is the learning rate. By default, γ(1) = 1. The parameters used by Shin & Chang for forgetting and learning were: γ = 0.0004, µ = 0.0001 η and β are estimated by:

ψ 0 (η(i)) -log η(i) = S 2 (i) -log S 1 (i) (8) β(i) = η(i) S 1 (i) (9)
The left part of the first expression is monotonically increasing, so it is simple to do reverse lookup in a table to find η 1

Implementation

This was the first method I tried since it is used in the source papers.

I implemented this solely in Matlab. I could not make the algorithm converge properly.

The main fault was that certain parameters would change only a very small amount over the course of the iterations, and in any case even when initialized to the correct value would still exhibit the same "drifting" behavior. I didn't pursue this method for more than a week before moving to other methods. Take 0.2 SISE Equations Method

Introduction

In an attempt to circumvent the problems encountered with the Gradient Descent method, I sought other works in the literature which treated the problem of estimating the parameters. I found a paper published by Song in 2008 which purported to show a new method which was fast and globally convergent. This method used special equations called SISE to estimate the parameters using the Newton-Raphson wero-finding algorithm.

Theoretical Overview

(Re)-Introduction of the Generalized Γ Distribution and the traditional approach to its solution

The Generalized Γ Distribution pdf used in this method differs from the one used in previous methods:

PDF of Generalized Γ Distribution:

f (x; θ) = βx βλ-1 σ βλ Γ(λ) exp -( x σ ) β . forx∈R + := (0, ∞) (10)
ML estimator of GΓD: θn := (σ n , βn , λn ) satisfies the following equations:

σn := n i=1 X βn i n λn 1 / βn λn := β n n i=1 X βn i log(X i ) n i=1 X βn i -1 βn n n i=1 log X i -log n i=1 X βn i + log(n λn ) -ψ 0 ( λ) = 0
Where σn is the shape parameter, We see that σ n is an expression in terms of β n and λ n ; λ n is expressed in terms of β n Solve the single equation resulting from substituting λn in (2) into (3) (3) is ill-behaved and difficult to solve numerically

SISE Equations

We will call the true parameter vector θ 0 = (σ 0 , β 0 , λ 0 ) Probability & expectation with respect to f (x; θ): P 0 and E 0 .

Let M 0 : t → E 0 {X t } First derivative: Ṁ0 Define: Ṙ0 . = Ṁ0 M0 and Ṙ0 (0) . = lim t→0 Ṙ0 (t) From this and a few other equations not shown here, we develop:

S(β) . = log(M 0 (2β)) -2 log(M 0 (β)) -log(1 + β( Ṙ0 (β) -Ṙ0 (0))) = 0 (11) T (β) . = M 0 (2β) M 2 0 (β) -(1 + ( Ṙ0 (β) -Ṙ0 (0))) = 0 (12)
Since we generally do not know β 0 , we can use instead a sample mapping:

M n : t → n i=1 X t i /n (13)
For the derivative we have: Ṁn (t) = n i=1 X t i log(X i )/n which estimates: Ṁ0 (t) = E 0 {X t } log X and:

Ṙn . = Ṁn (t) M n (t) (14) Ṙn (0) . = lim t→0 Ṙn (t) = n i=1 log(X i )/n (15)
the sample SISE equations are:

S(β) . = log(M n (2β)) -2 log(M n (β)) -log(1 + β( Ṙn (β) -Ṙn (0))) = 0 (16) T (β) . = M n (2β) M 2 n (β) -(1 + ( Ṙn (β) -Ṙn (0))) = 0 (17)
Estimation of Parameters with unknown shape index parameter λ 0

Step 1: Compute the estimator βn by means of NR method based on either ( 16) or (17)

Step 2: Compute the estimators σn and λn by using βn found in the previous step using the expressions σn = M n ( βn ) λn and λn = ( βn ( Ṙn ( βn ) -Ṙn (0))) -1 .

In regards to step 1, the NR root-finding equations are:

βn,k+1 = βn,k - S n ( βn,k ) Ṡn ( βn,k ) (18) βn,k+1 = βn,k - T n ( βn,k ) Ṫn ( βn,k ) (19)
Where

S n (t) = log(M n (2t)) -log(M n (t)) -log(1 + β( Ṙn (t) -Ṙn (0))) = 0; (20)
In the Appendix the author gives the derivative for S 0 :

Ṡ0 (t) = 2( Ṙ0 (2t) -Ṙ0 (t))- (21) Ṙ0 (t) -Ṙ0 (0) + t(R (2) 0 (t) -Ṙ2 0 (t)) t( Ṙ(t) -Ṙ(0)) + 1 (22)
I will simply assume we can use the same expression for n rather than 0.

Ṡn (t) = 2( Ṙn (2t) -Ṙn (t))- Ṙn (t) -Ṙn (0) + t(R (2) n (t) -Ṙ2 n (t)) β( Ṙ(t) -Ṙ(0)) + 1 (23)
S n (t) and Ṡn (t) depend on four terms:

Ṙn (t) = Ṁn (t) Ṁn (24) Ṙn (0) = lim t→0 Ṙn (t) Ṙ(2) n (t) = M (2) n (t) Ṁn Second derivative (25) Ṙ2 n (t) = Ṁn (t) Ṁn 2 (26) Square (27) F n (γ) = (28) ψ 0   1 N i = 1 N |x i | γ 1 N 2 N i=1 N j=1 |x i | γ log |xi| γ |xj | γ   + log 1 N 2 N i=1 j = 1 N |x i | γ log |x i | γ |x j | γ - 1 N N i=1 log |x i | γ = 0 (29)
We use the Newton-Raphson root-finding algorithm, which is:

γ(k + 1) = γ(k) - F n (γ(k)) Ḟn (γ(k)) (30)
taking care to reject negative values of γ and to monitor its progress so that it reduces F n . Shin & Chang note that the solution to (36) tends to diverge and cannot be obtained with approximately N < 400.

η and β are estimated by:

ψ 0 (η(k)) -log η(k) = S 2 (k) -log S 1 (k) (31) β(k) = η(k) S 1 (k) ( 32
)
where the statistics in question are given by the following:

S 1 = (1/N ) N i=1 |x i | γ ( 33
)
S 2 = (1/N ) N i=1 log |x i | γ (34)
The left part of ( 38) is monotonically increasing, so it is simple to do a reverse lookup in a table to find η.

Implementation

I implemented this method in Matlab. In spite of the apparent promise of the method as perceived from the paper, I could not cause this method to correctly converge. While convergence was acheived, the algorithm (apparently) did not converge to the correct parameters of synthetic data generated from gaussian and gamma distributions. However, in light of the fact that the pdf of this method has a different form than that of the others, and also considering the fact that perhaps at this point I had not learned to properly generate arbitrary pdfs using the inverse-CDF method, it is possible that I was mistaken in my assessment of this algorithm's accuracy. Thus this algorithm may merit a second, more rigorous evaluation, especially if it is seen to be faster and less prone to error than the methods more recently tried. Regardless of the success or failure of the Song method, it inspired me to try the Newton-Raphson method in a different way. I recalled that in the papers of Almpanidis and Shin & Chang the stated goal of the gradient descent algorithm was to find the zeros of the following equation:

F n (γ) = (35) ψ 0   1 N i = 1 N |x i | γ 1 N 2 N i=1 N j=1 |x i | γ log |xi| γ |xj | γ   + log 1 N 2 N i=1 j = 1 N |x i | γ log |x i | γ |x j | γ - 1 N N i=1 log |x i | γ = 0 (36)
We use the Newton-Raphson root-finding algorithm, which is:

γ(k + 1) = γ(k) - F n (γ(k)) Ḟn (γ(k)) (37)
taking care to reject negative values of γ and to monitor its progress so that it reduces F n . Shin & Chang note that the solution to (36) tends to diverge and cannot be obtained with approximately N < 400.

η and β are estimated by:

ψ 0 (η(k)) -log η(k) = S 2 (k) -log S 1 (k) (38) β(k) = η(k) S 1 (k) ( 39
)
where the statistics in question are given by the following:

S 1 = (1/N ) N i=1 |x i | γ ( 40
)
S 2 = (1/N ) N i=1 log |x i | γ (41)
The left part of ( 38) is monotonically increasing, so it is simple to do a reverse lookup in a table to find η.

Implementation

I implemented this method in Matlab. This method gave good results on synthetic data geerated using the inverse-CDF method in finding the correct values of the parameters. However it has the disadvantage of being so slow as to be unusable.

Newton-Hessian Method 0.4.1 Introduction

In talking with Nizar Bouguila about the previous method, he recommended that I use the Newton method with the Hessian matrix to solve for the three parameters η β, and γ simultaneously. The first step in implementing this was to derive the second derivatives of the log-likelihood function. Dr. Bouguila and I did this independently and he later confirmed my results.

Theoretical Overview

As a reminder, the log-likelihood for the generalized gamma distribution is:

L = log f x (x; η, β, γ) = log γβ η 2Γ(η) + (ηγ -1) N i=1 log |x i | -β N i=1 |x i | γ
We want to maximize this expression through a numerical solution.

Newton-Raphson method with Hessian matrix: Let θ = (η, β, γ). Apply the following recurrence relation:

θ k+1 = θ k - ∂L ∂θ k ∂ 2 L ∂ 2 θ -1
Gradient: Vector of first derivatives with respect to each parameter in θ. Hessian: Matrix of second derivatives with respect to each parameter in θ.

∂L ∂η = N (log β -ψ 0 (η)) + γ N i=1 log |x i | ∂L ∂β = N η β - N i=1 |x i | γ ∂L ∂γ = N N γ + η N i=1 log |x i | -β N i=1 |x i | γ log |x i | ∂ 2 L ∂η 2 = -N ψ (1) (η) ∂ 2 L ∂β 2 = -N η β 2 ∂ 2 L ∂η∂β = N β ∂ 2 L ∂β∂γ = - N i=1 |x i | N log |x i | ∂ 2 L ∂η∂γ = N i=1 log |x i | ∂ 2 L ∂γ 2 = -N γ 2 -β |x i | γ (log |x i |) 2 (42)

Implementation & Practical efforts

The basic method was implemented in Matlab. This method worked quite well on synthetic data generated using the inverse-CDF method. It was relatively fast and accurate. However, numerical problems were encountered when applying the method to real speech data feature vectors and even when applying it to concatenated series of synthetic feature vectors generated by inverse-CDF. The latter is explained as follows:

Generation of Synthetic Feature Data

Our speech feature vectors were of dimension 10 (10 DCT coefficients.) We had analysis window shifts of 10 ms and statistical windows of 500ms to 2s, giving us 50 to 200 feature vectors fo each statistical window and half that for each subwindow. For an experiment, we generated 30 GΓ Distributions. Starting with three base distributions -Gaussian, Gamma, and Laplacian, we generated 10 distributions from each base distribuion by perturbing the parameters with a random uniform number of appropriate scale. Then we generate a vector of random numbers drawn from each of the 30 distributions. This vector is of length approximately 3 times the length of a subwindow, or approximately 750ms to 3 s ==¿ 75 features to 300 features. There are ten such vectors for each base distribution, so we end up with three matrices of size 10x75 to 10x300, each drawn from random perturbations of the base distributions. In other words, we have one matrix corresponding to 10 variations of the Gaussian, another with 10 variations of the Gamma, and another with 10 variations of the Laplacian.

Based on these synthetic data, we perform the fitting of 10 GΓDs and find the log-likelihoods and GLRT distances.

Numerical Problems

It was in doing this that we ran into numerical problems. Namely, the local estimate of the parameters would occasionally tend to having two of the parameters of very large scale (¿50) with the other parameter approaching zero. This caused the Hessian matrix to become singular and the log-likelihood to become infinite in scale (probably an artifact of the hugely inflated parameters.)

Another numerical problem encountered was that, in the course of the iterations of the Newton algorithm, sometimes taking the full Newton step led to a negative η β, or γ. The log likelihood function cannot handle negative data or negative parameters.

Another problem is that sometimes, especially with small sample sizes encountered when using small subwindows, the algorithm was ill-behaved, wherein taking the full Newton step did not increase the log-likelihood but rather decreased it.

These three problemes were tackled using an implementation of backtracking as I understand it from the book "Numerical Recipes".

Outline of Newton-Hessian Algorithm as Implemented in Matlab

The algorithm is written as a Matlab function. It accepts as input the data and three control parameters. The first control parameter, λ, controls the scale of a term used to determine whether to run backtracking in a given iteration. The second, ρ, is used to control the scale of a measure which determines whether to stop the algorithm. The third, d, was once used as a control to stop the algorithm.

The Newton method is run within a while loop. The conditions for exiting the loop are the following:

(1)

2 |diff • diff| < 1 × 10 -6 or (2) # iterations > maxits
The program starts with parameters η β, and γ set to 1, 1, and 1, respectively. The algorithm takes the full Newton step computed as α ∂ 2 L ∂ 2 θ -1 ∂L ∂θ k , where by default α is set to 1 at the beginning of each loop/iteration. Next, the program goes into a while loop which checks if any of the parameters has gone negative with the step taken. If so, it records the current value of α as α 0 , and applies α = α -α 0 /10. It then re-estimates the parameters and if any is still negative, repeats until either the parameters are all positive or it reaches a maximum-loop-iteration counter. This limit is usually set at 25.

The next backtracking routine is as follows.

1. The program computes the following: term = λ ∂L ∂θ k (θ k+1 -θ k ) Originally, I had set this term as ∂L ∂θ0 , where this is the initial gradient on the first step of the algorithm. However this did not seem to work as well so I changed it to the current form.

2. α 0 = α 3. The program enters into a while loop. The loop runs while the updated log-likelihood J k+1 is less than or equal to J k -term, or while J k+1 is infinite or NaN. The loop terminates when these conditions are no longer satsified or when it reaches a limit on a loop counter. While the conditions of the loop are active, it repeats much the same procedure as the backtracking routine used to prevent the parameters from becoming negative: α = α -α 0 /10. Then we compute temporary parameters according to:

θ temp = θ k -α ∂ 2 L ∂ 2 θ -1 ∂L
∂θ k and compute a temporary updated loglikelihood J temp . If J temp is in fact less than the previous estimated log likelihood taking the full Newton step (i.e. J k+1 ,) or if J temp is infinite or NaN, or if any temporary parameter is less than zero, we break out of the while loop and discard the temporary parameters and temporary updated log-likelihood. Otherwise, J k+1 = J temp and θ k+1 = θ temp (the updated log-likelihood and parameters get overrwitten with the new values taking a smaller step.)

We repeat until J k+1 > J k -term or J k+1 is not infinite or NaN.

Previously, these backtracking routines did not subtract a fixed of α, but rather divided α by 2. This, however, could not resolve some cases. One backtracking routine which used this method and is no longer used in the present form of the NR-Hessian algorithm worked as follows:

1. We check if the updated likelihood J k+1 is less than J k and if, at the same time, J k <= J k-1 .

2. If so, we apply α = α/2.

θ

temp = θ k-1 -α ∂ 2 L ∂ 2 θ -1 ∂L ∂θ k-1
4. Re-check the first condition and if true repeat steps 2 and 3.

After the loop completes, we apply the following 1. θ k = θ temp .

Note that in this step we step back one step and revise θ k by taking a fractional Newton step from θ k-1 .

2. We revise J k accordingly 3. Erase θ k+1 and J k+1 4. Decrement k by 1.

After much use and modification, and more importantly after modification to the other backtracking routines, I felt this routine was unneccessary and it has been commented-out in recent revisions of the NR-Hessian algorithm.

Conclusion

I have still not been able to force this algorithm to work perfectly with the speech data. Nizar Bouguila suggested the use of the gammaln Matlab function rather than taking the logarithm of gamma(η) in the log-likelihood expression. This was supposed to help avoid infinite log likelihoods, but it doesn't seem to have helped much. The best I can do to prevent singular Hessian matrices and parameters and log-likelihoods which go to infinity is to emplace a condition which tells the while loop to quit if the magnitude of the Hessian becomes too great. However any such limit may prove too high depending on the kind of data presented. :

An example of the kind of data with which this algorithm has trouble. The blue bars represent the empirical pdf. The colored lines are successive fits of the GΓD. The black line is the final fit.

An example of the kind of data with which this algorithm excels Numerical problems aside, the algorithm performs well on synthetic data drawn from idealized GΓD distributions when the window size is rather large (not on the order of what would used for phoneme segmentation.) they are non-negative and returns the log-likelihood of those "fudged" parameters.

The second problem is not easily tractable. We can try to avoid it by limiting the number of function evaluations and by raising the tolerance which ends the loop. This will, however, theoretically result in a loss of accuracy.

The large-scale method performs approximately 15% slower than the mediumscale method.

Results with the Matlab Toolbox Methods

Synthetic Data

We tested the algorithms on the same synthetic data setup used before. We started with large window sizes applied to unidimensional data.

GLRT distance computed from log-likelihoods of GΓD's fitted by the fminunc method.

Statistical window size is 275 features.

Next we moved to smaller statistical window size and smaller segment lengths.

GLRT distance computed from log-likelihoods of GΓD's fitted by the fminunc method to speech data. The blue curve is the speech signal. The cyan lines are the manually transcribed phoneme boundaries. The green line is the GLRT distance. The red line is the same distance passed through a low-pass filter.

Log likelihoods of the fitted GΓD's for the large window and the two subwindows for speech data.

0. 3

 3 Newton-Raphson Algorithm for Parameter Estimation 0.3.1 Introduction & Description of Method

GLRT distance computed from log-likelihoods of GΓD's fitted by the fminunc method.Statistical window size is 60 features.We also applied a low-pass filter for this step.Low-pass filtered GLRT distance. Statistical window size is 60 features. Segment length is 360 features.

GLRT distance computed from log-likelihoods of GΓD's fitted by the fminunc method to synthetic data.Log likelihoods of the fitted GΓD's for the large window and the two subwindows.

Matlab Optimization Toolbox Methods

Several Matlab Optimization Toolbox methods were tried. These are divided into variations on the use of fminunc and the use of fminsearch. fminsearch is a non-derivative based secant method, whereas fminunc uses either a usersupplied analytical or finite-difference numerical gradient with the option of a user-supplied analytical or finite-difference numerical Hessian. fminunc is divided into two modes of operation: Medium-Scale and Large Scale. Medium-Scale uses a Quasi-Newton method, whereas Large Scale uses a Trust-Region Newton method.

Perfomance

While these methods all perform equivalently in terms of accuracy on synthetic data generated from idealized GΓD's with a large number of samples, we once again encounter numerical problems with the Medium-Scale fminunc and with fminsearch. We encounter two very familiar scenarios:

1. the algorithm leads the parameters to be negative 2. the algorithm leads the parameters to be very large resulting in infinite Log-Likelihood

The first problem is largely resolved through a fix in the routine which estimates the Log-likelihood. This adds a small increment to the parameters until Low-pass filtered version of GLRT distance for 10-dimensional feature synthetic data.

The increase in the number of dimensions seems to increase the signal-tonoise ratio of the distance graph.

The configuration of settings we settled on is as follows. We use the Large-Scale algorithm and supply the gradient as a function, but not the Hessian matrix. We reduce the maxiumum number of function evaluations from 300 to 150, increse the Function Tolerance ("Function" here meaning both the objective function of the log-likelihood and the first-order optimality) to 2e-3, and increase the tolerance on X (the parameters of the GΓD being estimated) to 5e-3. We take these last three steps in an effort to reduce the time taken to fit a distribution to the data.

Speech Data

I have applied the method to speech data to obtain the GLRT distances. On one sample of speech, the result is the following distance graph.


Caroline Ngo

Damien Koenig

Olivier Sename

Hubert Béchart

A Reduced Model of Three Ways Catalyst Converter and Stored Oxygen Rate Estimation Using Switched Observer

Pollutant emissions limitations ruled by the latest standards force car manufacturers to improve after treatment systems. This improvement is achieved by new systems on one hand, and by more efficient control strategies, on the other. This paper describes a 0-D model of Three Ways Catalyst converter (TWC) which is able to reproduce the transient dynamics of the converter although simple. A switched observer of stored oxygen rate is derived from this reduced model, and its convergence is proved. An estimation of the oxygen capacity of TWC provides precious information for on-board diagnosis or monitoring purpose.

First, a description of the TWC dynamics and its reduced model will be described, then a discrete time switched observer of the stored oxygen rate and the free site capacity estimation by least squares method are presented.

I. INTRODUCTION

Incomplete fuel combustion in cylinder leads to pollutants emissions. Nowadays, severe norms impose car manufacturers strict pollutant emissions regulations. As a result, exhaust gas after treatment systems need to become more efficient. Manufacturers meet the constraints by developing new technologies, but also through more sophisticated engine control strategies. The present study deals with the major after treatment system of Spark Ignition (SI) engine : the Three Ways Catalyst (TWC) converter.

The distinctive feature of Three Ways Catalyst is its ability to convert the three major regulated pollutants (CO, HC and N O x ). Pollutants optimum conversion require the fuel air ratio to be at stoichiometry value. Above this value, mixture is defined as rich (excess of fuel in mixture) whereas below, the mixture is defined as lean (excess of oxygen). The fuel air ratio, through the injected fuel quantity in the engine, is controlled by Exhaust Gas Oxygen (EGO) sensor which indicates the presence of oxygen in exhaust pipe. However, catalyst system control strategy based on exhaust fuel air ratio measurements, without taking account of its internal dynamics, is not enough to meet futures restrictions. Over the years, researches focused on TWC converter mathematical models. However, modeling this system is a complex task due to the numerous chemical reactions involved, as presented in [START_REF] Montreuil | Modeling current generation catalytic converters : Laboratory experiments and kinetic parameter optimization -steady state kinetics[END_REF] [START_REF] Pattas | Transient modeling of 3-way catalytic converters[END_REF].

Caroline Ngo and Hubert Béchart are with Renault, 1 Avenue du Golf, TCR LAB 0 12, 78288 Guyancourt Cedex, France (caroline.ngo, hubert.bechart)@renault.com Damien Koenig and Olivier Sename are with GIPSA-Lab, Grenoble INP -Université Joseph Fourier, 11 rue des Mathématiques, BP46, 38402 SAINT MARTIN D'HERES Cedex , France (damien.koenig, olivier.sename)@gipsa-lab.grenoble-inp.fr Many reduced models of TWC converter have been developed to satisfy new engine management strategy ([3][4][5] [11][17]). Reduced model from [START_REF] Brandt | Dynamic modeling of a three-way catalyst for si engine exhaust emission control[END_REF] is focused on the oxygen storage and release rates, described by a static function. Transport time delay was assumed to be negligible. This model is very simple, however, dynamic response to TWC transients is not taken into account.

As in [START_REF] Brandt | Dynamic modeling of a three-way catalyst for si engine exhaust emission control[END_REF], authors in [START_REF] James | A simplified model for the dynamics of a three-way catalytic converter[END_REF] described a simplified model dominated by the oxygen storage and release dynamics modeled by a polynomial function which parameters are identified from experimental data obtained at different engine conditions. The presented model include the effects of space velocity. These models are focused on the main dynamic of TWC : the oxygen storage dynamics. They are simple but have many parameters to adjust. This paper aims to respond to model based needs to meet future standards and present an application of switched observer to chemical process.

A 0-D physical model based on simplified dynamics of TWC from [START_REF] Aimard | Mathematical modeling od automotive three way catalytic converters with oxygen storage capacity[END_REF] is presented, also focused on oxygen storage dynamics. This model has been validated on experimental results. A switched linear observer has been designed for oxygen storage rate monitoring, designed by solving multiple Linear Matrix Inequality (LMI) and the obtained gains guarantee the stability whatever switching sequence. Estimation of oxygen free sites density is done using least squares algorithm to take into account this parameter deviation in observer. An application will illustrate the presented method.

Notations : (.) T stands for transpose matrix and (.) > 0 denotes a positive definite matrix.

II. MAIN DYNAMICS OF THREE WAYS CATALYST CONVERTER

A. TWC dynamics focused on oxygen storage

Exhaust gas of a SI engine is mainly composed of unburnt hydrocarbons (HC), carbon monoxides (CO), nitrogen oxides (N O x ), which are the major regulated pollutants, and oxygen (O 2 ), essential to convert pollutants. The TWC is able to convert these three pollutants when gas flows through it. Indeed, thanks to its internal structure covered by precious metal (Palladium, Platinum, Rhodium), the TWC promote HC/CO oxidation and N O x reduction reactions and convert them into water, carbon oxide and nitrogen.

N O x reduction or HC/CO oxidation are realized according to the exhaust gas composition characterized by the Fuel Air Ratio, denoted R. When the gas is rich (R > 1), N O x are eliminated by HC and CO present in excess, whereas, during lean phase (R < 1), HC and CO are converted by N O x present in excess this time. As illustrated by TWC efficiency curves (Fig. 1), best conversion rate for both reduction and oxidation reactions is located at stoichiometry value (R = 1).

FAR / FARstoich

Conversion rate (%) Fig. 1. Conversion efficiency However, keeping R at stoichiometry value is impossible due to the high transients in driving conditions. To increase the efficiency of TWC conversions, Ceria oxides (Ce 2 O 3 ) are added to fix oxygen. During lean phase, O 2 and N O x excess facilitate Ceria oxidation, once stored, they can be released to convert pollutants when mixture is rich. By varying the fuel air ratio around stoichiometry, oxygen storage and release are stimulated and pollutants conversion is possible for any air-fuel mixture if storage management is optimum. Efficiency is then between 90% and 100%, and also depends on TWC temperature and gas flow rate. Indeed, heat activate the catalyst (chemical kinetics) and flow rate, if too high, prevent gas to stay long enough for being convert (mixture transport along converter).

TWC is affected by sulfur dioxides poisoning, but also by heat damages (overheat during misfires) which decrease the efficiency of conversions in time. As a result, the density of free sites for oxygen storage decrease, less oxygen is stored and pollutants at exhaust increase. Ageing of TWC represented by density of free sites has to be taken into account in the model and estimated.

B. Chemical reactions

As described in the previous section, five reactions mainly describe TWC dynamics. Exhaust gas elements react with each other according those main chemical reactions :

-

N O x reduction : 1 x N O x + CO → CO 2 + 1 2x N 2 -HC oxidation : HC + 1 2 O 2 → 1 2 H 2 O + CO 2 -CO oxidation : CO + 1 2 O 2 → CO 2
The oxygen storage dynamics is modeled by the following reactions :

-Oxide adsorption when R < 1 :

1 2 O 2 + Ce 2 O 3 → 2CeO 2 -Oxide desorption and CO conversion when R ≥ 1 : CO + 2CeO 2 → CO 2 + Ce 2 O 3

C. Reduced model of Three Ways Catalyst

Previous reactions can be grouped to reduce the model : one involve reducers and the other, involve oxides. The exhaust fuel air ratio is expressed with concentrations of equivalent reducer and oxide by :

R out = 1 + K( 1 2 C out RED -C out OX ) (1)
All these simplifications lead to the following reaction rate equations :

-Oxygen rate dynamic : Other assumptions can be considered : when mixture is rich, it can be supposed that few oxides are present and adsorption is negligible, so k ads = 0 ; on the contrary, when mixture is lean, k red = 0 since reduction reactions are negligible.

θ = 2k ads C out OX (1 -θ) -k red C out RED θ -Exhaust fuel air ratio dynamic : Ṙout = vvh(R in -R out ) + Kρ 2 2k ads C out OX (1 -θ) -k red C out RED θ With : θ : Oxygen ratio stored in TWC R in ,
Expressions of oxide and reducer concentration from (1) can be reduced again and approximated by :

C out OX = 1-Rout K , if R out < 1, 0 otherwise C out RED = 2(Rout-1) K , if R out ≥ 1, 0 otherwise
The model is composed of 3 inputs

(u = [R in , vvh, T ]), 2 states (x = [θ, R out ]) and measured output (y = R out ).
The non linear system is given by :

R out < 1    ẋ1 = 2 K k ads (u 3 ) (1 -x 1 )(1 -x 2 ) ẋ2 = u 2 (u 1 -x 2 ) +k ads (u 3 ) ρ (1 -x 1 )(1 -x 2 ) R out ≥ 1    ẋ1 = 2 K k red (u 3 ) x 1 (1 -x 2 ) ẋ2 = u 2 (u 1 -x 2 ) +k red (u 3 ) ρ x 1 (1 -x 2 )
(2) Simulation of the non linear model results against experimental data is presented in the following figures : Results of the non linear model simulation shows good fit of the model to data (Fig. 2). The reduced model behave as expected : during rich phase, all stored oxygen have been released and θ = 0, whereas during lean phase, θ = 1 meaning that all free oxygen sites are occupied. The oxygen storage and release dynamics (TWC breathing) are also well modeled : during rich to lean and lean to rich transitions, R out = 1. In the first case, TWC store all the oxygen present in gas and since mixture is lean, no oxide and reducer are present are TWC exhaust, so R out = 1 during storage. In the second case, stored oxygen is released and react with reducer, as a result R out = 1 during this process.

In order to apply a switched observer, this reduced model is linearized.

III. SWITCHED LINEAR MODEL

The choice of switched system obtained by linearization over other form of representation is due to the nature of TWC's operating range (rich or lean ranges). In addition, the model is clearly non linear but dynamics are simple.

Let consider the following switched system :

   ẋ(t) = n i=1 µ i (x 2 (t))(A i x(t) + B i u(t) + d i ) y(t) = Cx(t) (3)
With :

A i = ∂f (x,u) ∂x xi,ui , B i = ∂f (x,u) ∂u xi,ui , C i = C = [0 1], d i = f (x i , u i ) -A i x i -B i u i And : ẋ(t) = f (x(t), u(t))
The linearized system is given by the following matrices (A i , B i , d i ) :

-If R out < 1 :

A i = 2 K k adsi (x 2i -1) 2 K k adsi (x 1i -1) ρ k adsi (x 2i -1) -u 2i + ρ k adsi (x 1i -1) B i = 0 0 2 K E ads ru 2 3i k adsi (x 1i -1)(x 2i -1) u 2i u 1i -x 2i E ads ru 2 3i ρ k adsi (x 1i -1)(x 2i -1)
d i = f (x i , u i ) -A i x 1i x 2i -B i   u 1i u 2i u 3i   -If R out ≥ 1 : A i = 2 K k redi (1 -x 2i ) -2 K k redi x 1i ρ k redi (1-x 2i ) -u 2i -ρ k redi x 1i B i = 0 0 2 K E red r u 2 3i k redi x 1i (1-x 2i ) u 2i u 1i -x 2i E red r u 2 3i ρ k redi x 1i (1-x 2i ) d i = f (x i , u i ) -A i x 1i x 2i -B i   u 1i u 2i u 3i  
Each model i represents an operating point and to switch between the n models, measured x 2 = y = R out is chosen as switching variable. µ i is the switching signal and satisfies the following convex properties :

   n i=1 µ i (x(t)) = 1 0 ≤ µ i (x(t)) ≤ 1
Here, µ i is constant and equals to 0 or 1 according to switching variable value so, if

µ i (x 2 (t)) = 1, (A i , B i , d i ) are activated.
Four set of points are chosen to model the dynamics (n = 4), they represent lean and rich steady state phases and fuel air ratio transitions. Operating points are obtained by mathematical identification based on non linear model data.

The chosen points are given in the following table :

x 1i x 2i u 1i u 2i u 3i Rout < 1 i = 1 0.

TABLE I SETS OF OPERATING POINTS FOR THE FOUR SUB MODELS

In the case where switching variable is not measurable, this switched linear model can be extended to estimate the switching variable. This approach is presented in [START_REF] Ichalal | State estimation of takagi-sugeno systems with unmeasurable premise variables[END_REF], [START_REF] Ichalal | Observer design and fault tolerant control of takagi-sugeno nonlinear systems with unmeasurable premise variables[END_REF]. Good results are obtained by this switched model although simple : modelization error of FAR remains small in spite of noise from switches. Whereas, model of oxygen rate from switched model is equivalent to the one from non linear model, and is not affected by switching sequence. Noise from switches of FAR is not a problem since this value is measurable. This model has been applied on data different from the one used in identification process which justify the validity of this switched model.

In the following sections, a discrete time switched observer and a parameter estimation are presented for x 1 = θ and ρ monitoring. Parameter estimation is used as an input for the observer in order to take into account the TWC ageing.

A. Switched observer synthesis

The expression of the discrete time switched observer is of the following form :

   xk+1 = n=4 i=1 µ i (x 2k )(A di xk + B di u k + d di + L i (y k -ŷk )) ŷk = C d xk (4
) Where L i is the observer gain of the current model i and A di , B di , d di and C d are the discrete matrix of A i , B i , d i and C, respectively (with a sample time of 10 ms).

These matrices depend on the known parameter ρ estimated by least squares algorithm described in Section V.

From (3) and ( 4), the state estimation error dynamic is given by :

e k+1 = n=4 i=1 µ i (x 2k )(A i -L i C)e k (5)
e k = x k -xk
Problem : Considering the switched model (3) and the switched observer (4), find the gains L i such that the estimation error dynamics ( 5) is globally asymptotically stable.

Theorem 1 : Observer gain For i, j ∈ 1, 2, ..., n, where i represent k and j, k+1, if the pair (A di , C) is detectable for all i, if there exist matrices P i ∈ R m×m , positive definite and symmetric and, U i ∈ R m×p such that the following LMI is verified for all (i, j) :

-P i A T di P T i -C T d U T i P i A di -U i C d P j -2P i < 0
then, the discrete time switched observer (4) of system (3) exist and the gains are given by L i = U i P -1 i .

Sketch of proof : (steps similar to [START_REF] Daafouz | Stability analysis and control synthesis for switched systems : a switched lyapunov function approach[END_REF], [START_REF] Koenig | Unknown input observers for switched nonlinear discrete time descriptor systems[END_REF])

To guarantee convergence toward 0 of the estimation error (5), the following inequality should be satisfied :

V (k + 1, e k+1 ) -V (k, e k ) < 0
for all k, where V (k, e k ) = e T k P k e k > 0 is a candidate Lyapunov function, with P k = P T k and P k > 0.

Let P i = P k and P j = P k+1 , the difference V (k + 1, e k+1 ) -V (k, e k ) along the solution of (5) becomes :

-P i +(A T di P T i -C T d L T i P T i )P -1 i P j P -1 i (P i A di -P i L i C d ) < 0
By Schur's complement, this inequality is equivalent to :

-P i (A di -L i C d ) T P T i P i (A di -L i C d ) -P i P -1 j P i < 0
With U i = P i L i , and by noticing that P T i P -1 j P i ≥ 2P i -P j , the LMI to solve becomes :

-P i A T di P T i -C T d U T i P i A di -U i C d P j -2P i < 0

B. Switched observer synthesis with pole assignment

A solution to improve the performances of the observer is to add constraints on matrices A di and C d in the previous LMI to force poles assignment in sub-region ( [START_REF] Koenig | Unknown input observers for switched nonlinear discrete time descriptor systems[END_REF], [START_REF] Vinicius | Discrete-time switched systems : Pole location and structural constrained control[END_REF]). r , where σ is the small circle center coordinates and r is the small circle radius.

The LMI to solve becomes :

-P i ( A di -σI r ) T P T i -( C d r ) T U T i P i ( A id -σI r ) -U i ( C d r ) P j -2P i < 0

C. Observer results

The switched observer has been validated against the same data used previously in switched model validation paragraph. Initial states values are different from data to see observer performances. Pole placement application will not be illustrated.

Figures show that observer converge and is not affected by switches. Maximum value of error estimation for x1 and x2 is 30%. Although this value is quite high, error converge quickly to 0 such as the results are satisfactory.

The presented results prove that observer is efficient for oxygen storage rate monitoring (Fig. 7) which can be employed in advance control strategies including oxygen storage management.

V. FREE SITES DENSITY ρ ESTIMATION BY LEAST SQUARES ALGORITHM

Future standards will impose On Board Diagnosis (OBD) of TWC. Catalytic converters are mostly affected by sulfur poisoning or overheat which reduce the surface available to convert pollutants. As a result, density of free oxygen site ρ decrease, leading to efficiency loss of TWC. By monitoring the parameter ρ, the TWC ageing can be taken into account by the observer to adapt oxygen storage management, and is a precious information for OBD.

To monitoring ρ, an estimation by least squares method is chosen (see [START_REF] Landau | System Identification and Control[END_REF]). One can remark that an adaptive observer, PI observer or addition of the parameter in the state vector could have been used to estimate both states and parameter as in [START_REF] Guyader | Adaptive observer for discrete time linear time varying systems[END_REF] or [START_REF] Chadli | State and unknown input estimation for discrete time multiple model[END_REF].

However, two reasons explain this choice :

-The parameter ρ varies slowly (years), so online estimation is not useful for observer. -System with addition of ρ as a state and ρ(k+1) = ρ(k), is not observable when x 1i = θ i = (0, 1) or x 2i = 1 (eq.( 2)). As a result, decoupling the problem of observation and estimation was preferable.

The recursive least squares estimation is done with the values from the non linear model. The non linear model presented in section II is of the form :

h(x) = ρf (x, u) + g(x, u) With : h(x) = ẋ2 g(x, u) = u 2 (u 1 -x 2 ) f (x, u) = k ads (u 3 )(1 -x 1 )(1 -x 2 ) (if x 2 < 1) f (x, u) = k red (u 3 )x 1 (1 -x 2 ) (if x 2 ≥ 1)
Since ẋ2 is used to estimate the parameter, the least square algorithm is more advantageous : data from non linear model are stored to run the algorithm afterward. Moreover, the parameter ρ varies slowly, the sample time of estimation process can be smaller than the observer one (see Fig. 5) to provide known parameter to observer as online estimator.

The objective is to estimate the parameter ρ using least squares criteria : The following result is obtained by running estimation algorithm with a sample time of Ts = 1e-4 s with an observer running with a sample time of Te = 0.06 s. The estimation converge to the parameter value even when parameter real value changes, showing that this algorithm used in a recursive way is fast and for monitoring application, especially the ageing of TWC.

J Θ (N ) =

VI. CONCLUSION

In this paper, a reduced model well suited for control thanks to its simplicity in the form of ordinary differential equations of small order, is developed based on oxygen dynamics of the TWC.

A multiple model linearization representation permits to design a discrete switched observer by solving multiple LMIs, considering all the possible cases. Structural constraints can be added in observer synthesis to force its poles to be in a specified circle. A parameter estimation based on recursive least squares algorithm is provided, treated separately from the observer problem.

Successful application of the oxygen fraction rate θ observation and estimation of the free site density ρ will improve the oxygen storage management, needed for future standards.

  HC/CO/H 2 can be represented by an equivalent reducer, denoted RED, and O 2 /N O x by equivalent oxide, denoted OX. Indeed, in first approximation, CO and H 2 have similar behavior and quantities of HC and N O x are negligible compare to other components of exhaust gas. The previous chemical equations become : -Oxygen adsorption : OX + 2Σ → 2Σ * -Pollutants reduction : Σ * + RED → Σ + CO 2 /H 2 O With : Σ/Σ * : Free/occupied oxygen storage site OX : Equivalent oxide RED : Equivalent reducer

Fig. 2 .

 2 Fig. 2. Oxygen rate model and Fuel air ratio from experimental measurements (dashed) compared to the non linear model (solid)

  Validation is done from non linear model data obtained by simulation with experimental data in input. Inlet FAR is a step signal between the values of R in = 0.92 and R in = 1.05 with non constant sample time, temperature and mass flow rate are constant. Model identification results are shown in the following figures :

Fig. 3 .

 3 Fig. 3. Fuel air ratio and Oxygen fraction rate from non linear model (solid) compared to the switched model (dashed)

Fig. 5 .

 5 Fig. 5. Observation and estimation process scheme

1 Fig. 6 .

 16 Fig. 6. Sub-region for pole location

Fig. 7 .

 7 Fig. 7. Fuel air ratio and Oxygen fraction rate from non linear model (solid) compared to observer result (dashed)

  i) -Θ T Φ(i -1)) 2With : y : Measured output Φ : Vector composed of measured output and input Θ : Vector of parameter to be estimated Recursive least squares algorithm can be applied by writing :y(x, u) = h(x) -g(x, u) Φ(x, u) = f (x, u) Θ = ρ

Fig. 8 .

 8 Fig. 8. Free site density estimation (solid) -Parameter value (dashed)