Keywords: Adaptive Innovation-based method, Fault diagnosis, Adaptive amplifier coefficient, Automotive bicycle model

This article is concerned with a fault detection enhancement method using adaptive amplifier coefficients (a.a.c.) concept. It is developed for linear systems with mixed uncertainties (stochastic and bounded uncertainties framework). The study provides, also, analysis and discussions about the applicability and the efficiency of the enhanced method to several sensors fault error types. Simulations on a vehicle bicycle model (validated by experimental tests on the Renault Megane) are presented to emphasize on the performances of the developed method.

Fault detection using Enhanced Adaptive degrees of freedom χ 2 -statistics method for Linear systems with mixed uncertainties Quoc Hung Lu * Soheib Fergani * Carine Jauberthie *

INTRODUCTION

Within the control theory and its field of applications, fault detection is extremely important for all system engineering problems. It is also a crucial component of any system diagnosis scheme and has received a lot of attention recently. Reliable fault detection and isolation (FDI) is a first class requirement in many fields. Indeed, efficient (early and accurate) fault detection can help avoid dangerous scenarios (accidents, explosions,...) or improve productivity (reducing process activity loss such as leakage...). The model-based approaches are proven to provide good results and acceptable tradeoff between fault sensitivity and computational cost especially those based on residual generation (see [START_REF] Patton | Issues of fault diagnosis for dynamic systems[END_REF] and references within).

Many methods and techniques have been developed to meet these abundant requirements. Several methods for fault detection in dynamic systems are mentioned in [START_REF] Willsky | A survey of design methods for failure detection in dynamic systems[END_REF], including the innovation-based method in which a χ 2 -statistic hypothesis testing was used. This method is applied appropriately with the standard Kalman filter [START_REF] Kalman | A new approach to linear filtering and prediction problems[END_REF] to process the linear dynamic system with deterministic coefficient matrices. In [START_REF] Sainz | Fault detection and isolation of the three-tank system using the modal interval analysis[END_REF], an approach to generate envelopes based on interval techniques of the modal interval analysis is proposed. [START_REF] Puig | Fault diagnosis and fault tolerant control using set-membership approaches: Application to real case studies[END_REF] reviews the use of set-membership methods in fault diagnosis and fault tolerant control. It aims at checking the consistency between observed and predicted behaviour by using simple sets to approximate the exact set of possible behaviours. Also, the design of stable interval observers for linear systems with additive time-varying zonotopic input bounds is proposed in [START_REF] Raka | Fault detection based on robust adaptive thresholds: A dynamic interval approach[END_REF]. Interval state observers provide an estimate on the set of admissible values of the state vector at each instant of time. Ideally, the size of the evaluated set is proportional to the model uncertainty, thus interval observers generate the state estimates with estimation error bounds, similarly to Kalman filters, but in the deterministic framework. Main tools and techniques for design of interval observers are reviewed in this tutorial for continuous-time, discrete-time and time-delayed systems (see [START_REF] Efimov | Design of interval observers for uncertain dynamical systems[END_REF] )). In [START_REF] Efimov | Design of interval observers for uncertain dynamical systems[END_REF] a survey on interval state observers that provide an estimate on the set of admissible values of the state vector at each instant of time is provided. The efficiency of these strategies has attracted the attention of the industrial community. The automotive industry, especially, has given more attention to these methods and though many academic studies have tried to provide solutions within this field based on set membership FDI. In [START_REF] Meseguer | Fault diagnosis using a timed discrete-event approach based on interval observers: Application to sewer networks[END_REF] a fault diagnosis approach is proposed. It has been motivated by the problem of detecting and isolating faults of the Barcelona's urban sewer system limnimeters (level meter sensors). It is based on interval observers which improve the integration of fault detection and isolation tasks. [START_REF] Ifqir | Robust interval observer for switched systems with unknown inputs: Application to vehicle dynamics estimation[END_REF] reviews the problem of robust state estimation and unknown input interval reconstruction for uncertain switched linear systems. A design method for obtaining interval observers that provide guaranteed lower and upper bounds of the state and unknown inputs is applied to vehicle lateral dynamics estimation to show the effectiveness of the algorithms. Also, in [START_REF] Chen | Robust estimation of vehicle motion states utilizing an extended set-membership filter[END_REF], an extended set-membership filter applied to the vehicle's longitudinal velocity, lateral velocity, and sideslip angle provides not only higher accuracy, but also can provide a 100% hard boundary which contains the real values of the vehicle states (compared to the Unscented Kalman Filter UKF-based approaches). Recently, in [START_REF] Tran | Cadre unifié pour la modélisation des incertitudes statistiques et bornées -Application à la détection et isolation de défauts dans les systèmes dynamiques incertains par estimation[END_REF], authors provided an approach based on a χ 2 -statistics test whose degrees of freedom (d.f.) is fixed and predetermined thanks to a size of a widows in which the statistic is computed. This approach is an innovation-based method but it is modified and applied to interval residual context. Then, [START_REF] Lu | A new scheme for fault detection based on optimal upper bounded interval Kalman filter[END_REF] have proposed an adaptive degrees of freedom χ 2 -statistic (ADFC) method to provide an improved solution for the same problem (mixed uncertainty assumption, interval residual context) of the former. This development consists in investigations of interval analysis related to d.f. of the considered statistic and results in the use of the so called adaptive amplifier coefficients (a.a.c.).

In present paper, the main contribution is the development of an enhanced ADFC (eADFC) method by discussing different levels of a.a.c. used in the method. The Bicycle vehicle model, used as a benchmark, is presented in detail together with multiple simulation results of fault detection based on this model.

The paper is organized as follows. Section 2 provides useful notations, the linear discrete time dynamic system with potential sensor faults, the enhanced method also evaluation indicators. In Section 3, the application to fault detection based on the Bicycle vehicle model is presented, included detail presentations of the model, simulation settings and results. Section 4 is a conclusion.

ENHANCED METHOD

Notations

A real interval matrix [X] of dimension p × q is a matrix with real interval components [x ij ], i ∈ {1, ..., p}, j ∈ {1, ..., q}. Write X ∈ [X] to indicate a point matrix X = (x ij ) belonging element-wise to [X]. Define: X ≡ sup([X]) △ = (sup([x ij ])) , X ≡ inf([X]) △ = (inf([x ij ])), as element-wise operators applying to [X] and then mid([X]) △ = (X + X)/2, rad([X]) △ = (X -X)/2, width([X]) △ = X -X. Define also the (convex) hull of two interval matrices [X 1 ], [X 2 ] of the same dimension as hull{[X 1 ], [X 2 ]} △ = [min{X 1 , X 2 }, max{X 1 , X 2 }].
Basic interval operators ⋄ ∈ {+, -, ×, ÷} defined in [START_REF] Jaulin | Applied Interval Analysis, with Examples in Parameter and State Estimation, Robust Control and Robotics[END_REF] can be used to compute directly all operations [u] ⋄ [v] and α ⋄ [u], for real intervals [u], [v] and α ∈ R, without any further approximation algorithm. Then, interval matrix computations are defined similarly to matrix computations using the basic operators while more general operators are constructed by means of inclusion function [f ] [START_REF] Jaulin | Applied Interval Analysis, with Examples in Parameter and State Estimation, Robust Control and Robotics[END_REF]. In practice, the package Intlab [START_REF] Rump | Intlab -interval laboratory[END_REF] developed for Matlab (also existing in Octave and C/C++) is used for computations.

Denote I(x) = 1 if the conditions x are true and I(x) = 0 otherwise, mean(x) = n i=1 x i /n, ∀x ∈ R n , and S + ([M ]) the set of all positive semi-definite matrices belonging to an interval matrix [M ]. The notation p : l : q (p ≤ q) is used for a range from p to q with step l. For l = 1, we write p : q. A sequence of variables can be noted interchangeably as w 1 , ..., w k or w 1 : w k or w 1:k .

Model

Consider the following linear discrete time dynamic system

x k = A k x k-1 + B k u k + w k , y k = C k x k + D k u k + v k + f s k , k ∈ N * , (1)
where x k ∈ R nx and y k ∈ R ny represent respectively state variable and measure,

u k ∈ R nu input, w k ∈ R nx and v k ∈ R ny stochastic noises. f s k ∈ R ny is a sensor fault vector.
Each of its components corresponds to a sensor fault. Thus, the fault vector f s k can be of the multiple or single error type. In the first type, some (or all) sensors cause errors which affect the y k value for the corresponding components. In the second type, only one sensor causes errors and just the corresponding y k component is affected.

Assumptions (H):

The initial state x 0 is normally distributed with mean µ 0 and covariance P 0 . w k and v k are centered Gaussian vectors with covariance and[µ 0 ] respectively. x 0 , {w 1:k } and{v 1

Q k and R k . Matrices A k , B k , C k , D k , Q k , R k , P 0 and µ 0 are unknown and belong to given interval matrices [A], [B], [C], [D], [Q], [R], [P 0 ]
:k } are mutually independent.
System (1) with assumptions (H) is a quite general model adapted to a wide range of applications. In this system, parameter matrices are time varying, the uncertainty may result from different sources (system disturbances, measurement noises) and may be of different kinds (stochastic and bounded uncertainties).

To evaluate the fault detection performance, indicators introduced in Lu et al. ( 2021) are used. Assume that system (1) is applied for N iterations among which faults occur in an time interval R with length l (0 ≤ l ≤ N ). R may be an interval or union of them but is called hereafter an error range for simplicity. A detection signal is denoted by π k . It has value 0 (no fault is detected) or 1 (a fault is detected). A right detected signal is a signal π k = 1 with k ∈ R and a false detected signal is a signal π k = 1 with k ̸ ∈ R. Then, the indicators are defined as follows:

+ Detection Rate: DR = k∈R I(π k = 1)/l × 100%, + False Alarm Rate: FAR = k̸ ∈R I(π k =1)
N -l × 100%, + Efficiency: EFF = DR -FAR, + No detection rate: NDR = 100% -DR.

Enhanced Adaptive degrees of freedom χ 2 -statistic (eADFC) method

The ADFC in Lu et al. ( 2021) is an innovation-based method but deals with interval residual

[r k ] = y k -[C][x k|k-1 ] -[D]u k . [r k ]
is the difference between the obtained measure y k and an interval estimate of the later. [r k ] contains all admissible residuals r k 's computed according to assumptions (H).

[x k|k-1 ] = [A][x k-1|k-1 ] + [B]u k is the interval version of the a priori estimate xk|k-1 of the (standard) Kalman filter. The considered statistic is U k (Algorithm 1) which can be expressed as a function of width([r k ]). Furthermore, U k is an upper bound of all admissible χ 2 statistics r T k S -1 k r k of n y d.f
. where S k is the covariance of r k . Therefore, [START_REF] Lu | A new scheme for fault detection based on optimal upper bounded interval Kalman filter[END_REF] considered U k as a χ 2 statistic of κ k n y d.f. with an a.a.c. chosen as

κ k = mean {width([r k ])} .
(2) This enlarges the traditional χ 2 statistic whose d.f. is normally a positive integer. The enlargement is mathematically legal thanks to the continuity of the χ 2 cumulative distribution function in its argument and parameter.

In general, we claim that :

• the choice of κ k is not unique, • there exists a better/optimal choice of κ k for each of different applications (due to its intrinsic conditions),

• any choice of κ k depends on predetermined purposes, e.g. control the system using u k as control input or maintain the FAR below a target (5% for instance).

Therefore, in the present paper, we do not aim to provide a generic solution of the optimal choice of κ k , nor to solve this problem in some specific cases. Instead, a flexible enhanced method is proposed by controlling the κ k strength using a level parameter λ k > 0. The controlled a.a.c. is now expressed as:

κ k = λ k • mean {width([r k ])} .
(3) The method is simple and meaningful, since the κ k of (2) is shown in [START_REF] Lu | A new scheme for fault detection based on optimal upper bounded interval Kalman filter[END_REF] to be chosen reasonably according to analysis of the statistic U k , i.e. κ k is sensitive to the fault, large enough to obtain a small FAR (e.g. ≤ 5%) in the fault free case, sensitively affected by width([r k ]) as well as U k but it does not increase as fast as the later whenever a fault occurs and affects on the width([r k ]). This choice also provides a good performance fault detection with FAR < 3% in several scenarios of simulations. Thus, the κ k of ( 3) is raised naturally and its different controlled levels can be easily adapted to different applications and purposes.

In addition, we let open the choice of a.a.c. for designers but provide some guidelines and propositions as follows: P1) Constant level. This is the most simple use of the level parameter, in which λ k ≡ c, ∀k ≥ 1, and c > 0 is determined experimentally by simulating many scenarios of the considered application. The best value of c according to some criteria is chosen. P2) Piece-wise constant level. This is an advanced version of P1, in which λ k ≡ c i , for N i < k ≤ N i + M , c i > 0, {M, N i , i} ⊂ N and N 0 = 0. The levels c i can be determined such that FAR computed after every M time instants is bounded in some range [x, x]. Once FAR is smaller than x the level must be increased and vice-versa. Other criteria can be proposed to determined the level parameters.

The approach P1 is necessarily implemented off-line in the design stage while the approach P2 is apt to be implemented online to provide updated λ k 's after every delay of M time instants. This online procedure can be effectuated by using simultaneously two eADFC detectors as shown in Fig. 1. Indeed, the first detector aims at generating the detection signals for the system. The second one consists in providing an updated level λ k by computing FAR indexes and using virtual faults f test Finally, the adjusting procedure of Algorithm 1 is important to enhance the fault detection method. It can help to reduce the FAR index by considering that in a window of size W (before or after time instant k), if the number of consecutive detected faults is not beyond W , then these faults (if exist) don't cause serious effects and will be dismissed. More precisely, π k is kept its value 1 if k+j i=k-W +1+j π i = W for j = 0 or j = W -1, which is expressed equivalently by line 14 of Algorithm 1.

The enhanced method is summarized and presented by Algorithm 1 in the following.

Algorithm 1 eADFC method

1: Initialization: [x 0|0 ], P 0|0 , [A], [B], [C], [D], [Q], [R], α, λ k , W , u k , y k , k = 1, 2, ..., N . 2: for k = 1, 2, 3, ...N do 3:
Compute: [x k|k-1 ] and [P k|k-1 ] using OUBIKF (Algo.1 in [START_REF] Lu | A new scheme for fault detection based on optimal upper bounded interval Kalman filter[END_REF]).

4:

[r k ] = y k -[C][x k|k-1 ] -[D]u k 5: [S k ] = [C][P k|k-1 ][C] T + [R] 6:
U k = sup{abs([r k ] T [r k ]/a k )} 8: κ k = λ k • mean {width([r k ])} 9: Find δ k s.t.: P(χ 2 (κ k n y ) > δ k ) = α. 10: Detection signal : π k = I(U k > δ k ). 11: end for Adjusting procedure: 12: Let π k=1:N +2W ← [0 1×W , π k=1:N , 0 1×W ] 13: for k = 1 + W, ..., N + W do 14: π k ← or I k i=k-W +1 π i = W , I k+W -1 i=k π i = W
15: end for 16: π k=1:N ← π k=1+W +⌊W/2⌋:N +W +⌊W/2⌋ Notes: 0 < α < 1, λ k > 0, W ≥ 1.

APPLICATION

In this section, the efficiency of the proposed method is highlighted using a vehicle model validated on a real car.

Bicycle model

The vehicle model parameters obtained by an identification process on the Renault Mégane Coupé are presented. Indexes i = {f, r} and j = {l, r} are used to identify vehicle front, rear and left, right positions, respectively. The full vehicle model with the nonlinear equations describing its dynamical behaviour can be found in [START_REF] Fergani | Robust multivariable control for vehicle dynamics[END_REF].

Since the full model is highly nonlinear, a linear bicycle model as illustrated by Fig. 2 reproducing the lateral behaviour of the car is used for this study by linearizing the former. Reference to Fig. 2, β(t) is the sideslip angle and ψ(t) is the vehicle yaw which form the model state variables. F ty f (t) represents lateral front tire forces, F ty r (t) represents lateral rear tire forces and F tx f (t) represents the longitudinal front tire forces, v is the vehicle speed, ∆F txr(t) is the differential rear braking force (obtained based on the braking torques T brj ), δ is the steering angle and M dz is the yaw moment disturbance.

The model is obtained considering the following: a) Low sideslip angles: |β| < 7 degrees, b) Low longitudinal slip ratio: < 0.1, c) Low steering angles: cos(δ) ≃ 1. The linearized lateral tire forces are:

F ty f (t) = C f β f (t), F ty r (t) = C r β r (t), (4)
with the front and rear sideslip angles β f (t) and β r (t):

β f (t) = δ(t) -β(t) - l f ψ(t) v , β r (t) = β(t) + l f ψ(t) v . (5)
This leads to the following state space representation (6):

β(t) ψ(t) =    -C f -C r mv 1 + µ -l r C r -l f C f mv 2 -l r C r -l f C f I z -l 2 f C f -l 2 r C r I z v    β(t) ψ(t) +    C f mv 0 0 0 l f C f I z 1 I z S r Rt r 2I z - S r Rt r 2I z       δ M dz T b rl T brr    .
(6) Remark 1. The sideslip dynamics are highly nonlinear. They cannot be measured thanks to conventional sensors. µ is the tire/road adhesion coefficient. It takes value in [0, 1] depending on the road conditions (dry, wet, icy,...) and influences strongly the vehicle lateral dynamics.

Simulation

The discretization of the model is performed with a time period T = 0. Inputs u k 's are simulated according to a dynamic change for N = 864 time instances, i.e. the vehicle is assumed to be driven at 15 m/s on a dry road (µ = 1) and a double line change maneuver is performed from t = 0.5s to t = 1.5s by the driver. The initial state is chosen at x 0 = (0, 0) T . For every k, choose at random A comparison between the proposed method and a fixed d.f. χ 2 fault detection strategy shows the importance of the adaptive characteristic of the approach. The method proposed in Tran (2017) (called Method A in the following) uses the statistic

A k , B k , C k , D k , Q k , R k in corresponding interval matrices and so that Q k ⪰ 0 and R k ⪰ 0. Then, w k ∼ N (0, Q k ), v k ∼ N (0, R k
T k = inf([r k ] T S -1 k [r k ]) with S + ([S k ]) ⪯ S k .
In this method, a fault is detected if T k > δ where δ is determined by P(χ 2 (W n y ) > δ) = α. The first disadvantage of method A is that interval computation can let T k be negative, consequently no fault is detected as illustrated in Fig. 4 according to the bicycle model simulation. The second disadvantage is that a windows size W is arbitrarily chosen. An example can be built to illustrate this method works quite well (Fig. 5 Left) in which T k is non negative, but then the second disadvantage is still critical: another choice of W leads to another detection result. Our adaptive method in this case still provides an accurate fault detection (Fig. 5 Right). For this sake, the following scenario setting are used. Let f s k 's be identical error vectors with k in an error range R whose length is l = 50. For each of different values of f s k , choose randomly error range R and do L = 100 times of y k simulations and fault detections. Indicators are computed for each of L simulation times and their means are yielded afterward. Let τ k = max{f s k }/Max width where max{f s k } is the maximum component of f s k and Max width is the maximum width of diagonal elements of [Q] and [R]. This quantity provides an idea of how large is the maximum fault value with respect to some known quantity propagating through the dynamic system and affecting the measure values y k , that is the maximum variance of noises. Remark 2. This scenario is a combination of those scenarios presented in [START_REF] Lu | A new scheme for fault detection based on optimal upper bounded interval Kalman filter[END_REF], in which the error range or the y k simulation is fixed. The following remarks are valid for all three types of errors already simulated:

(R1) FAR does not vanish even in the fault free case (b = 0). This fact implies that there are other reasons (than fault) causing FAR. Actually, in this case, the error range degenerates to length 0, all Tables 4 and5 show that FAR increases when the level parameter λ k decreases from 0.7 to 0.3. EFF also rises up for small values of b (< 35) according to the decrease of c.

When c = 0.7 (Table 4), the values of the FAR index disperse in [3.8 , 5.1](%), and thus EFF = DR -FAR does not overpass 96.2% even if DR reaches its maximum value (100%). In addition, compared to Table 2 (c = 1), EFF increases considerably for many fault values b < 35 (starting at 10) and begins to achieve a remarkable rate (64.11%) starting at b = 25. Also in comparison with Table 2, but then the EFF does not increase in the case b = 5; this is due to the fact that Tables 2 and4 display the simulation results of the different samples. The case b = 0 is not comparable (see remark (R1)).

When c = 0.3 (Table 5), the FAR range is [23.9 , 25.3](%), EFF is never beyond 76.1%.

So, depending on the applications requiring a low FAR or a high EFF for a fine fault detection (detecting error with small fault value), different a.a.c κ k can be chosen suitably thanks to the level parameters λ k > 0.

CONCLUSION

An enhanced adaptive method to fault detection is proposed in this paper. The enhancement is presented in connection with the former method. The simplicity and efficiency of the method is highlighted from the design idea to the implementation perspective (as shown in the application simulations). The proposed method is well combined with the active fault diagnosis (AFD) approach developed in [START_REF] Lu | Active fault diagnosis based on adaptive degrees of freedom x2-statistic method[END_REF].

Fig. 1 .

 1 Fig. 1. General eADFC piece-wise constant level implemention diagram

Find

  a k s.t. : S + ([S k ]) ⪯ a k I using Theorem 2 in Lu et al. (2021).

Fig. 2 .

 2 Fig. 2. View of the bicycle model reproducing the lateral behaviour of the car.

  05s. Then, point matrices A, B, C, D are obtained and interval matrices [A], [B], [C], [D] are generated as follow: for M ∈ {A, B, C, D}, let M = mid([M ]) and choose the radii rad([M ]) at random in [0 , 0.5]. The covariance matrices [Q] and [R] are generated in the same way by noting that their diagonal elements must be intervals of positive real numbers.

  ) are simulated and {x 1:N , y 1:N } are computed according to (1). Sensor faults are generated in terms of additive vectors f s k ∈ R ny . Let b, b ′ ∈ R be constant fault values, following types of error will be handled in the remaining: Type 1: f s k = b.1 where 1 is the all-ones vector in R ny . Type 2: f s k = b.e j where e j is the j-th column of the identity matrix I ny . Type 3: f s k = b.e j + b ′ .e j ′ , with j ̸ = j ′ . The error terms are added to y k for all k in a range R with length l, i.e. R = r : r + l -1 for r ∈ {1 : Nl + 1}. Every sequence {y 1i , y 2i , ..., y N i }, with i ∈ {1, ..., n y }, is called a chain. So, the errors occurred on multiple chains of y k in type 1 and type 3 and only on single chain j in type 2. Algorithm 1 is applied with following initializations: α = 0.05, W = 5, [x 0 ] = ([-0.5, 0.5], [-0.5, 0.5]) T , P 0|0 = max{diag(sup([Q]))}I = 0.4412 I , upper bounds ω k I of any set S + ([M ]) identified by ω k = ∥Max([M ])∥ F (Frobenius norm) where [M ] = ([m ij ]) and Max([M ]) = (max ij ) so that max ij = sup([m ij ]) if mid([m ij ]) ≥ 0 and max ij = inf([m ij ]) otherwise.

3. 3

 3 Results using identity level parameter λ k ≡ 1 First results. A result of fault detection using eADFC method for type 1 of error with fault value b = 20 is shown in Fig. 3. The error range [722 : 772] is between the two vertical black lines. The detection signals are very well determined.

Fig. 3 .

 3 Fig. 3. eADFC method-Fault detection for Bicycle model.

Fig. 4 .

 4 Fig. 4. Method A -Fault detection using statistic T k = inf([r k ] T S -1 k [r k ]) for Bicycle model.

Fig. 5 .

 5 Fig. 5. Example-Fault detection with bias sensor fault b = 10: (Left)Method A, (Right) eADFC method. Advanced results. Now, in order to survey, using eADFC method, how well the detection is when influencing factors are changed (e.g. fault values b and b ′ , error range R and simulated variable y k ), indicators introduced in section 2 are used in following simulations. Three types of errors are simulated to get insight of the proposed method.

  Type 1 error. The sensor faults come to all chains of {y k }, f s k = b.1 and τ k = τ = b/Max width with b = 0 : 5 : 30.

  Type 2 error. The sensor faults only occur in one chain of {y k }, f s k = b.e j and τ k = τ = b/Max width with b = 0 : 5 : 60. The chain j on which the faults occur is chosen randomly at each of L times of y k simulations. This situation corresponds to a single sensor and normally not all sensors are damaged at the same time. This situation is necessary for fault isolation in a further phase.Type 3 error. Let b = 10.m and b ′ = 10.(m + 1) for m = 1, 2, 3. The chains j and j ′ at which the faults occur are also chosen randomly at each simulations. This setting, while still being of the multiple error type, can represent an intermediate situation between the settings of type 1 and type 2 error previously presented.

Table 1 .

 1 Fault detection for type 1 error.

	b	τ	DR%	NDR%	FAR% EFF%
	0	0	1.48	98.52	1.53	-0.05
	5	13.8	4.54	95.46	1.59	2.95
	10 27.5	13.88	86.12	1.69	12.19
	15 41.3	72.38	27.62	1.97	70.41
	20 55.0	95.38	4.62	1.75	93.63
	25 68.8	99.62	0.38	2.13	97.49
	30 82.5	99.96	0.04	1.96	98.00

Table 2 .

 2 Fault detection for type 2 error.

	b	τ	DR%	NDR%	FAR%	EFF%
	0	0	3.28	96.72	1.46	1.82
	5	13.8	3.74	96.26	2.01	1.73
	10	27.5	3.58	96.42	1.44	2.14
	15	41.3	10.26	89.74	1.66	8.60
	20	55.0	26.50	73.50	1.74	24.76
	25	68.8	42.26	57.74	1.88	40.38
	30	82.5	50.40	49.60	1.80	48.60
	35	96.3	75.58	24.42	1.64	73.94
	40 110.0	78.08	21.92	1.52	76.56
	45 123.8	85.22	14.78	1.72	83.50
	50 137.5	89.88	10.12	2.91	86.98
	55 151.3	96.24	3.76	2.10	94.14
	60 165.0	98.18	1.82	2.05	96.13

Table 3 .

 3 Fault detection for type 3 error.

	(b, b ′ )	DR%	NDR%	FAR% EFF%
	(10,20)	42.84	57.16	5.29	37.55
	(20,30)	79.38	20.62	5.33	74.05
	(30,40)	98.82	1.18	5.59	93.23

Table 4 .

 4 In the following simulations, λ k ≡ c, ∀k ≥ 1, with c ∈ {0.7, 0.3}. The simulation results are shown in Tables 4-5 in which only the type 2 error is concerned. Both cases of c are applied for common data samples. Fault detection for type 2 error using a.a.c κ k with level parameter λ k ≡ 0.7

	However, we can think that b has a very small value
	and thus results in the Tables remain unchanged.
	(R2) The negative value for EFF at b = 0 in Table 1 can
	be explained by the fact that, in this case, the fault
	detection procedure not only provides no efficiency
	gain, but rather a loss.
	(R3) The factors causing FAR are multiple. Two of these
	factors that differ from simulations to simulations
	are random noises and random error ranges, which
	can therefore be called specific factors. Some other
	factors that exist for all simulations, and which can
	therefore be called general factors, are: the model
	performance (how well the model describes the dy-
	namics of the vehicle), the interval computations, the
	lack of knowledge on the exact values of A k , B k , C k ,... and the performance of the χ 2 -statistic test
	(with α significance level).
	(R4) DR and EFF have ascending trends according to b
	values (or to max{b, b ′ }) while FAR is rather stable
	in some range with positive values.
	(R5) There is a threshold for good/bad result of EFF, e.g.
	b ≥ 15 for type 1, b ≥ 35 in type 2, max{b, b ′ } ≥ 30
	in type 3.
	(R6) The choice of significance level α is also a tunning
	factor for an appropriate fault detection.
	3.4 Results for other constant level parameters λ

1-value detection signals are false detected signals, DR and NDR are not defined and FAR must be recomputed, e.g. according the first row of Table

1

: FAR = [1.53 × (N -50) + 1.48 × 50] /N ≈ 1.527. k

In this section, some other constant level parameters applied to Bicycle model are investigated. So, the performance of the eADFC method with constant level (P1) can be completely illustrated. The piece-wise constant level case (P2) is built straightforwardly from the former.

Table 5 .

 5 Fault detection for type 2 error using a.a.c κ k with level parameter λ k ≡ 0.3

	b	τ	DR%	NDR%	FAR% EFF%
	0	0	24.64	75.36	25.14	-0.50
	5	13.8	37.12	62.88	24.52	12.60
	10	27.5	61.84	38.16	25.12	36.72
	15	41.6	84.50	15.50	24.08	60.42
	20	55.0	95.46	4.54	24.72	70.74
	25	68.8	99.42	0.58	24.83	74.59
	30	82.5	99.82	0.18	24.06	75.76
	35	96.3	99.84	0.16	23.96	75.88
	40 110.0	99.92	0.08	25.25	74.67


domain: spi chim sde

text:
Yendoube Charles

Sano Moyeme

Stéphanie Betelu

Johan Bertrand

Stephane Gaboreau

Karine Gronen-Serrano

Ioannis Ignatiadis

à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Influence of polysulfide radicals on the mechanisms and intensity of corrosion of a carbon

for more than eight months, into a metric corrosion pilot, in which the steel was immersed into the cementitious FS10B material (water/solid ratio of 3). The system used a solid AgCl/Ag reference electrode that showed its robustness over the period and a Platinum counter-electrode. In complement, potentiometric measurements (pH, temperature, conductivity, redox potential, permittivity, dissolved O2), and continuous analyzes of the level of sulfides were simultaneously made. Particular attention was given to the visual analysis of the evolution of the grout color, which acquired a green tint in the first days of the experiment, due to the anoxic conditions prevailing into the hydrated FS10B material and the presence of polysulfide radical anions (Sn .-) coming from the blast furnace slag used for CEM III/A cement production.

Indeed, the loss of the green tint was observed after five months of experiment. It suggested Sn .-oxidation in agreement with FS10B partial dehydration and oxygen diffusion into the cement material. In addition to EIS, this specific behavior helped to understand, in particular, the reduction reaction involved during the API 5L X65 steel corrosion process and to choose the best kinetic model to apply during the first five month of experiment and beyong.

During the first five months of immersion into the cementitious FS10B material, the API 5L X65 steel was in an active state (2-5 µA/cm²) in the presence of Sn .-, despite the formation of a layer of corrosion products highlighted by EIS (Fig. 1A), and the corrosion current density was limited by the diffusion of Sn .-. Beyond five months, in the absence of Sn . , the steel become passive and the corrosion was limited by charge transfer of water molecules. These changes occurred in a few days and the corrosion current dropped down to about 0.02-0.07 µA/cm² (Fig. 1B).

Whatever the experimental conditions, the corrosion kinetics of API 5L X65 immersed in FS10B is dominated by the cathodic reduction reaction of water molecules. The

  latter is catalyzed by the presence of Sn .-. The use of the specific electrochemical methodology developed for corrosion sensing of API 5L X65 immersed in FS10B provides redundant information: it clearly identifies the changes at the API 5L X65 steel/FS10B interface as well as both the corrosion mechanisms and the kinetics.

Figure 1 .

 1 Figure 1. A) Evolution of impedance measurements in Nyquist mode derived from EIS measurements on API 5L X65 electrode immersed in a cement grout as a function of time. Evolution of API 5L X65-grout interface with time. B) Evolution of corrosion current densities (CCD) of an API 5L X65 carbon steel immersed in a cement grout as a function of time. LPR: CCD determined by Linear Polarization Resistance, EIS: CCD determined by EIS, CA: CCD determined by Chrono-Amperometry.

steel API 5L X65 in contact with cement grout in future nuclear waste disposal program Influence of polysulfide radicals on the mechanisms and intensity of corrosion of a carbon steel API 5L X65 in contact with cement grout in future nuclear waste disposal program

	Yendoube Charles SANO MOYEME 1 2 , Stéphanie BETELU 1 , Johan BERTRAND 3 ,
	Stéphane GABOREAU 1 , Karine GRONEN-SERRANO 2 , Ioannis IGNATIADIS 1 ,
	1 BRGM (French Geological Survey), Orléans, France, 2 LGC (University Paul
	Sabatier, UT3), Toulouse, France, 3 ANDRA, (French national radioactive waste
	management agency), Châtelet-Malabry, France
	Cigeo is the future French deep Nuclear Waste Disposal (NWD) project for long-lived
	highly radioactive wastes. The clayey Callovo-Oxfordian formation (COx) was chosen
	as natural barrier to isolate the waste from people and environment. The waste will be
	placed in API 5L X65 steel linings, separated from the COx by a cemento-bentonic
	material (FS10B at the time), whose formulation may still evolve, and which will, among

other things, keep the steel in a passive state thanks to its basic pH (>12). For safety reasons and reversibility of storage, continuous monitoring of the Corrosion Current Densities (CCD) or corrosion rate of the carbon steel in contact with this cementitious material is therefore necessary. That is the framework of this work.

A specific electrochemical methodology, consisting of i) stationary techniques (measurement of polarization resistance by Linear Polarization Resistance (LRP) and by Chrono-Amperometry (CA), Voltammetry around Open Circuit (VaOCP) by fitting the kinetic parameters), and ii) transitory techniques (measurement of impedance by Electrochemical Impedance Spectroscopy (EIS)), was carried out at regular intervals,


Hung Quoc

Soheib Lu

Carine Fergani

Jauberthie

Reinforced Likelihood Box Particle Filter

Keywords: State filtering and Estimation, Nonlinear System, Box Particle Filter, Interval Analysis

This paper is concerned with the development of a general scheme for box particle filtering. It is based on the likelihood computation, the most crucial step of the estimation strategy. The proposed filter takes advantages from strong aspects of various existing box particle filters and adds an interesting reinforced likelihood computation method that enhances the estimation results. An overview on Box Particle Filters and discussions from assumptions used in the literature to the filters performance evaluation approach are presented. Also, comparative study of the obtained results by performing several scenarios on the illustration example is provided to highlight the efficiency of the proposed estimation strategy.

I. INTRODUCTION

I N State Estimation or Filtering problems, when dealing with a linear Gaussian state-space model, analytical expressions computing the state estimates according to posterior distributions can be derived by the well known and widespread Standard Kalman Filter (SKF) [START_REF] Kalman | A New Approach to Linear Filtering and Prediction Problems[END_REF]. Many extension of SKF are then provided by numerous researches in different contexts [START_REF] Mohamed | Robust finite-horizon kalman filtering for uncertain discrete-time systems[END_REF]- [START_REF] Lu | Optimally bounded interval kalman filter[END_REF]. For nonlinear model without Gaussian measurement assumption, Particle Filters (PF) have been applied successfully to a variety of state estimation problems [START_REF] Gordon | Novel approach to nonlinear/non-gaussian bayesian state estimation[END_REF], [START_REF] Doucet | Sequential Monte Carlo Methods in Practice[END_REF]. The PF efficiency and accuracy depend mostly on the number of particles used in the estimation which may require a high computation time.

One of the famous extensions of PF to set membership approach is the Box Particle Filter (BPF) [START_REF] Abdallah | Box particle filtering for nonlinear state estimation using interval analysis[END_REF]. BPF handles box (interval vector of) states and bounded errors by using interval computation and constraint satisfaction techniques. This method has been shown to control quite efficiently the number of required particles, hence reducing the computational cost and providing good results in several experiments.Therefore, it has been applied in many applications, included aerial [START_REF] Freitas | Autonomous crowds tracking with box particle filtering and convolution particle filtering[END_REF] and ground vehicle [START_REF] Tehrani Niknejad | On-road multivehicle tracking using deformable object model and particle filter with improved likelihood estimation[END_REF] estimations.

Since then, numerous variants of BPF are developed to deal with measurements bounded uncertainty, measurements stochastic uncertainty or measurements mixed uncertainty. Various techniques and theories have been proposed to address the diversity of requirements in these contexts, e.g. weight updating using Bayesian filtering technique extending to box particle case [START_REF] Blesa | State estimation and fault detection using box particle filtering with stochastic measurements[END_REF] or belief function theory with different methods [START_REF] Nassreddine | State estimation using interval analysis and belief-function theory: Application to dynamic vehicle localization[END_REF], [START_REF] Tran | Evidential box particle filter using belief function theory[END_REF].

In the present work, regarding this large variety of BPF, a scheme is proposed to give a generalized description that highlights the specificity of this class of filters. An analysis of the likelihood computation methodology is investigated. Theory Date of first submission for review: April 19, 2022. Authors are affiliated at LAAS-CNRS, Université de Toulouse, CNRS, UPS, Toulouse, France. Emails: qhlu,sfergani,cjauber@laas.fr background is also provided to support the development of a novel filter benefiting advantages of existing methods via a numerous reinforcement techniques: score function, reduction percentage, exponential weighting, backward estimate.

The paper is organized as follows. The problem formulation is presented in Section II with discussions about assumptions used in the literature. Section III presents the general scheme of BPF. The likelihood computation methodology is analyzed in Section IV which also provides the notion of likelihood distinguishability. The proposed method is developed in Section V and simulations of the method are provided in Section VI. Section VII is the paper conclusion.

II. PROBLEM FORMULATION

A. Notations and definitions

A real interval matrix [X] of dimension p×q is a matrix with real interval components [x ij ], i ∈ {1, ..., p}, j ∈ {1, ..., q}. Write X ∈ [X] to indicate a point matrix X = (x ij ) belonging element-wise to [X]. Define: X ≡ sup([X]) = (sup([x ij ])) , X ≡ inf([X]) = (inf([x ij ])),
and then mid(

[X]) = (X + X)/2, rad([X]) = (X -X)/2, width([X]) = X -X. Denote also X = sup([X]), X = inf([X]), [X] = [X, X] and define the (convex) hull of two interval matrices [X 1 ], [X 2 ] of the same dimension as hull{[X 1 ], [X 2 ]} = [min{X 1 , X 2 }, max{X 1 , X 2 }].
Basic interval operators ∈ {+, -, ×, ÷} defined in [START_REF] Jaulin | Applied Interval Analysis, with Examples in Parameter and State Estimation, Robust Control and Robotics[END_REF] can be used to compute directly all operations [u] [v] and α [u], for real intervals [u], [v] and α ∈ R, without any further approximation algorithm. Then, interval matrix computations are defined similarly to matrix computations using the basic operators while more general operators are constructed by meant of inclusion function [START_REF] Jaulin | Applied Interval Analysis, with Examples in Parameter and State Estimation, Robust Control and Robotics[END_REF]. In practice, the package Intlab developed for Matlab is used for these computations.

B. Assumptions and discussions

Consider the following dynamical system:

(Σ) : x k = f (x k-1 , u k , w k ) , y k = h(x k , u k , v k ) , k ∈ N * , (1)
where x k ∈ R nx and y k ∈ R ny are respectively state and measurement output,

u k ∈ R nu input, w k ∈ R
Σ k ∈ R ny×ny . (D2) µ k , Σ k belong to known intervals [µ k ], [Σ k ]. (D3)
The measurements are point values y k .

Assumption (A) is used in [START_REF] Abdallah | Box particle filtering for nonlinear state estimation using interval analysis[END_REF], [START_REF] Blesa | State estimation and fault detection using box particle filtering with stochastic measurements[END_REF]- [START_REF] Tran | Evidential box particle filter using belief function theory[END_REF]. Assumptions (B) are under study in [START_REF] Abdallah | Box particle filtering for nonlinear state estimation using interval analysis[END_REF], [START_REF] Nassreddine | State estimation using interval analysis and belief-function theory: Application to dynamic vehicle localization[END_REF]. In [START_REF] Abdallah | Box particle filtering for nonlinear state estimation using interval analysis[END_REF], the BPF is introduced and becomes standard for many extensions or variants with essential steps: initialization, box propagation, contraction, likelihood (weight) computation, state estimation and resampling. In [START_REF] Nassreddine | State estimation using interval analysis and belief-function theory: Application to dynamic vehicle localization[END_REF], the Belief State Estimation algorithm is developed using the belief function theory. It may require some techniques for the construction and computation of masses, but after being normalized, these masses become likelihoods in the probability sense. Therefore, we also call likelihood computation as an essential step of this method.

Assumptions (C) are used in [START_REF] Blesa | State estimation and fault detection using box particle filtering with stochastic measurements[END_REF]. The method proposed therein includes a different approach to weight the box particles as well as a resampling procedure based on repartitioning the box enclosing the propagated states. There is no contraction step in this method.

Assumptions (D) are used in [START_REF] Tran | Evidential box particle filter using belief function theory[END_REF], in which (D1) is a special case of (C1) with a slight relaxation by adding bounded uncertainties to Gaussian parameters µ k and Σ k . In [START_REF] Tran | Evidential box particle filter using belief function theory[END_REF], the belief function theory is used with continuous mass functions to represent these kinds of uncertainties and to compute box particle likelihoods. The proposed approach therein leads to the so-called Evidential Box Particle Filter (EBPF) including all essential steps of the standard BPF.

Remark 1: (B3) is the implicit assumption deriving the consistency between the predicted measurement boxes

[h]([x i k ], [u k ]), i ∈ {1, .
.., M } (M the number of partitioned boxes), and the real measurement box [y k ]. This consistency is used in the contraction step and the likelihood computation by penalizing all particle boxes with which the intersections

[h]([x i k ], [u k ]) ∩ [y k ]
are empty. Remark 2: Assumptions (D3) and (C2) are coincided. They can be transformed into (B2) with a slight relaxation of (B3). That is, knowing the density of v k , we deduce its confidence intervals [v k ] with some significant level α and define

[y k ] = y k -[v k ].
Then (B3) is relaxed in the sense that the observed measurements [y k ] do not contain h(x k , u k ) with certainty but with only a probability (1 -α).

III. GENERAL SCHEME OF BOX PARTICLE FILTER

In general, although applying different background theories, the proposed methods in [START_REF] Abdallah | Box particle filtering for nonlinear state estimation using interval analysis[END_REF], [START_REF] Blesa | State estimation and fault detection using box particle filtering with stochastic measurements[END_REF]- [START_REF] Tran | Evidential box particle filter using belief function theory[END_REF] study State Estimation in a framework of stochastic uncertainties and/or bounded uncertainties with two main objectives :

Objective 1: Reduce the width of box particles to penalize the conservatism due to interval computations.

Objective 2: Quantify (compute) box particle likelihoods, as best as possible, to enhance the accuracy of the estimates.

The methods used in these references can be considered as variants of BPF and be summarized by Scheme 1 which is applied in a mostly similar manner across them.

Scheme 1 General Scheme of Box Particle Filtering

STEP 1 : Initialization. At a time instance k 0 ≥ 0, (re)partition the interval [x k0 ] or resample the set {[x j k0 ], w j } j=1:N k 0 into M disjoints equal-volume subboxes with the same weights:

{[x i k0 ], w i = 1/M } i=1:M . while a predetermined Condition C is still satisfied do STEP 2 : Propagation. Get a new set of box particles {[x i k0+1 ] = [f ]([x i k0 ], [u k0 ])} i=1,.
..,M estimating the box containing the real state x k0+1 = f (x k0 , u k0 ) with or without a contraction step. STEP 3 : Likelihood computation a) Compute (and normalize) the likelihoods of box particles {[x i k0+1 ]} i=1:M being the box containing the real state x k0+1 . This computation bases on the consistency between the estimated measurement [h]([x i k0+1 ], [u k0 ])'s and the obtained measurement [y k0+1 ] using different criteria and methods. By this step, the following set of box particles with updated weights is obtained : {[x i k0+1 ], w i } i=1:M . b) Some techniques can be applied at this step to get a more "efficient" set of box particles, e.g. discarding the boxes with small weights (smaller than some predetermined threshold) and with/without replicating the box associated with the greatest weight,... From this, the set of box particles becomes

{[x i k0+1 ], w i } i=1:N k 0 +1 , 1 ≤ N k0+1 ≤ M . STEP 4 : Estimation Interval estimate : [x k0+1 ] = N k 0 +1 i=1 w i .[x i k0+1 ] (2)
Point estimate :

x k0+1 = N k 0 +1 i=1 w i .mid([x i k0+1 ]) (3)
STEP 5 : k 0 = k 0 + 1 end while STEP 6 : Restarting at STEP 1.

Remark 3: In this Scheme, for a general presentation, the observed measurements are denoted as intervals since the point values are considered as special cases of intervals. N k0 in the initialization step takes value in {1, ..., M }. It is the number of box particles obtained at the end of the likelihood computation step at the previous time instant (k 0 -1). For k 0 = 0, the initialization concerns only the partition of [x 0 ] and not the resampling. The Condition C of the while loop in the Scheme may differ depending on the used method, e.g. in [START_REF] Abdallah | Box particle filtering for nonlinear state estimation using interval analysis[END_REF], this condition is N k0 ≥ N ef f where N ef f is a threshold calculated at each iteration of the loop.

IV. LIKELIHOOD COMPUTATION ANALYSIS A. Likelihood computation methodology

The diagram in Fig. 1 is provided to explain the methodology of LCMs used in Scheme 1.

Criteria to decide how a box particle be chosen as preferable than others Likelihood Computation Method (LCM)

Information issued from assumptions The assumptions of the system under consideration supply the information, or Info for short, needed to build the likelihood. The information can be, for instance:

• Info (a): The intersection between [y k ] and the box [h]([x i k ], [u k ]) containing the real value y k must be non empty, ([h]([x i k ], [u k ]) ∩ [y k ] = ∅). • Info (b): The distribution of v k and hence of r k = y k - h(x k , u k ) is Gaussian (for additive noise v k ).
The information can be directly an assumption or a deduction of the later. In bounded-error context, only Info (a) is treated [START_REF] Abdallah | Box particle filtering for nonlinear state estimation using interval analysis[END_REF] while in the mixed uncertainty case, both Info (a) and Info (b) are taken into account [START_REF] Tran | Evidential box particle filter using belief function theory[END_REF]. Criteria and methods are then chosen to exploit the information. Once a criterion is chosen, different methods can be used to calculate the likelihood. Inversely, a calculation method may correspond to many criteria. More precisely, a likelihood computation criterion is an application of the supplied information to establish an order of preferability of related events/objects that we want to compute their likelihoods, while a LCM is a method using mathematical formulae to represent and quantify the preferability by giving scores and normalizing them. Definition 1 formalizes these concepts.

Definition 1: Any set {s i > 0} i=1:n can be a set of scores. This set has an increasing (decreasing) order of preferability if the greater (smaller) score is preferable. The set

{L i } i=1:n is a set of likelihoods if L i ∈ [0, 1], i ∈ {1, ..., n} and n i=1 L i = 1. A normalizing function is a function Φ : R n → [0, 1] n
transforming scores to likelihoods so that a preferable score corresponds to the greater likelihood.

Let's illustrate the above ideas via Example 1.

Example 1: To exploit Info (a), one can use either:

• Criterion 1: The particle [x i k ] giving a "bigger" intersection determined by

[ẑ i k ] = [h]([x i k ], [u k ]) ∩ [y k ] must be preferable, • Criterion 2: The particle [x i k ] making [ŷ i k ] = [h]([x i k ], [u k ]
) "closer" to [y k ] must be preferable. How to represent "bigger" (size) or "closer" (closeness) notions and to calculate the corresponding likelihoods depend on the choice of LCMs. Criterion 1 is used in [START_REF] Abdallah | Box particle filtering for nonlinear state estimation using interval analysis[END_REF]. The corresponding LCM uses the volume Vol(.) determined by

Vol([x]) = nx j=1 width([x j ]), [x] = ([x 1 ], ..., [x nx ]) T , (4)
to represent the box size and computes the likelihoods L 1 = (L 1 1 , ..., L M 1 ) by :

s i 1 = Vol([z i k ]) Vol([ŷ i k ]) , i ∈ {1, ..., M } , L i 1 = s i 1 M i=1 s i 1 , (5)
where {s i 1 } are scores with increasing order of preferability. However, other methods can be used to calculate the likelihoods, e.g. with scores s j = (s 1 j , ..., s M j ), j ∈ {2, 3, 4}:

s i 2 = Vol([z i k ]), s i 3 = width([z i k ]) ∞ , s i 4 = width([z i k ]) 2 ,
in which s 3 and s 4 use distances between the two bounds of the intersection to represent its size. Criterion 2 can be applied with different LCMs using a distance between [ŷ i k ] and [y k ] to represent their closeness. Criterion 2 can also be used to exploit Info (b) as in [START_REF] Tran | Evidential box particle filter using belief function theory[END_REF] via the central tendency of the Gaussian vector

[r i k ] = y k -[ŷ i k ]
along with the belief function theory. A more detail analysis of the method used in [START_REF] Tran | Evidential box particle filter using belief function theory[END_REF] is found in section IV-A.

It is worth to note that, in some cases, it is difficult to distinguish clearly between criterion and LCM as illustrated by Fig. 1, e.g. in [START_REF] Blesa | State estimation and fault detection using box particle filtering with stochastic measurements[END_REF] with Interval Bayes filtering approach or in [START_REF] Nassreddine | State estimation using interval analysis and belief-function theory: Application to dynamic vehicle localization[END_REF] and [START_REF] Tran | Evidential box particle filter using belief function theory[END_REF] with the belief theory. The reason is that the criteria are implied under complexes theories.

B. Distinguishability of computed likelihoods

In order to deal with Objective 1, in the literature, contractors are usually applied based on the Constraint Satisfaction Problem technique. However, this is not the most crucial step of BPFs using Scheme 1, e.g. this step is skipped in [START_REF] Blesa | State estimation and fault detection using box particle filtering with stochastic measurements[END_REF]. Furthermore, partition a box into M disjoint equal-volume sub-boxes and then compute the expected interval by (2) also help to reduce the conservatism due to interval computations. The most crucial step that differs one method from another in this class of BPFs is the Likelihood computation focusing on Objective 2. This is thus the main discussion of this section. Let's begin by introducing the definition of distinguishability of likelihoods (Definition 2). Definition 2: Let {L i } i=1:n be a set of likelihoods, define the distinguishability of two likelihoods of the set by δ i,j ≡ δ j,i = |L i -L j |, 1 ≤ i = j ≤ n, and the total distinguishability of the set by δ T = n-1 i=1 n j=i+1 δ i,j ≡ i<j δ i,j . In general, the box likelihoods are computed at every time instance k. The more they can represent the ability of a box containing the real value, the better estimate is obtained by Estimation step. A weak distinguishability means that most of the computed likelihoods are quasi equal and hence not useful for distinguishing between box particles. In the remainder of this section, two representative groups of criteria and LCMs used in the literature are analyzed to show their major disadvantage which is the weak distinguishability.

Group I : Apply Criterion 1 with LCMs using the box volume for the box size representation.

This criterion is used implicitly in the Contraction step of all BPF algorithms including it and applied in [START_REF] Abdallah | Box particle filtering for nonlinear state estimation using interval analysis[END_REF] with the LCM L 1 defined by [START_REF] Gordon | Novel approach to nonlinear/non-gaussian bayesian state estimation[END_REF].

Consider real interval vector [y k ] = [y k , y k ] and real point vector δ such that 0 ≤ δ ≤ y k (element-wise), δ ∈ R p . Let T = diag{t 1 , ..., t p } be a diagonal matrix with diagonal entries {t r ≥ 0} r=1:p . Then, all boxes [ŷ i k ]'s with forms:

y k -T δ, y k + δ or [y k -δ, y k + T δ]
have the same likelihoods L i 1 = 1/ p r=1 (1 + t r ) (using ( 5)). In this case, these likelihoods have null distiguishability.

There are many other cases in which likelihoods are quasi equal and thus making the corresponding boxes [ŷ i k ] weakly distinguishable. For instance, M boxes [ŷ i k ]'s may have likelihoods L i = 1/M + δ i , δ i ∈ (-, ) with a small > 0. In this case, the benefit of the likelihood computation step could be insignificant.

Group II : Use Criterion 2 to exploit Info (b). The LCM used in [START_REF] Tran | Evidential box particle filter using belief function theory[END_REF] is investigated as the representative method of this group to deal with stochastic or mixed uncertainties and with additive Gaussian measurement noises. In this method, the innovation term

r k = y k -h(x k , u k ) is Gaussian with µ k ∈ [µ k ] and Σ k ∈ [Σ k ]. It belongs to some of intervals [r i k ] = y k -[ŷ i k ], i ∈ {1, .
.., M }. One defines:

[HV α ] = µ k -αDiag(Σ k ), µ k + αDiag(Σ k ) ,
where α ≥ 0, (.) is an element-wise operator and Diag(X) returns the diagonal of matrix X as a vector. Then, the belief bel(.) and plausibility pl(.) of [r i k ] are computed and considered as lower and upper bound of the probability of [r i k ] containing the real value r k as follows: bel

([r i k ]) = F ny+2 (α i bel ) and pl([r i k ]) = 1 -F ny+2 (α i pl ), α i bel = max{α : [HV α ] ⊆ [r i k ]}, α i pl = min{α : [HV α ] ∩ [r i k ] = ∅}, bel([r i k ]) ≤ Probability [r i k ] r k ≤ pl([r i k ])
, where F ny+2 is the cumulative distribution function of the χ 2 distribution with n y + 2 degrees of freedom. At this stage, Criterion 2 is applied based on the central tendency of the Gaussian vector r k : the more

[ŷ i k ] is close to y k (equivalently [r i k ] is close to [µ k ]
), the greater belief and plausibility the box particles [x i k ] (that yield [ŷ i k ] via the function h) attains. Finally, the likelihood of each particle [x i k ] is computed thanks to the Generalized Bayes theorem (GBT) and Pignistic transformation [START_REF] Tran | Evidential box particle filter using belief function theory[END_REF].

The weak distinguishability of the method is shown via the following critical point. All boxes [r i k ] intersecting [µ k ] have the plausibility 1. So, these boxes are not distinctive regarding their plausibilities. They are distinguished only by their beliefs, in which:

• For the boxes that intersect [µ k ] but do not contain it, their beliefs are 0. A zero information can be issued about these boxes. • For the ones containing [µ k ], their beliefs are characterized by the greatest focal element [HV α i bel ] they contain. The greater [HV α i bel ] a box can contain, the more belief it gets. It is quite similar to apply the rule: "the more [r i k ] is centralized (having a bigger intersection with [µ k ]) and has a bigger volume, the greater likelihood it gets". Other LCMs can be applied using that rule with a lightened calculation strategy and background theory.

Therefore, the likelihoods computed in the next step using GBT and Pignistic transformation are weakly distinguishable. The computation formulae are as follows [START_REF] Tran | Evidential box particle filter using belief function theory[END_REF]:

L i k = A⊂Ω,A =∅ m(A|y k ) |A| • 1([x i k ] ∈ A), ∀[x i K ] ∈ Ω, (6) m(A|y k ) = η [x i k ]∈A pl [r i k ] [x j k ] / ∈A 1 -pl [r j k ] , (7)
where A is a subset of

Ω = {[x i k ], i = 1 : M } with cardinality |A|, [r i k ] = y k -[h] [x i k ], [u k ] , 1(x) = 1 if
x holds true and vanishes otherwise and η = 1

-[x i k ]∈Ω 1 -pl [r i k ] . Example 2: Let Ω = {[x 1 k ], [x 2 k ], [x 3 k ] ∈ IR 2 } and put [r i k ] = y k -[h]([x i k ], [u k ]), i = 1, 2, 3. Assume that v k ∼ N (µ k , Σ k ), µ k ∈ [µ k ], Σ k ∈ [Σ k ], [r i k ] ∩ [µ k ] = ∅, i = 1, 3, and [r 2 k ] ∩ [µ k ] = ∅. Then, we get pl([r i k ]) = 1 for i = 1, 3 and 0 < pl([r 2 k ]) < 1.
Using ( 6) and ( 7), we obtain :

L 1 k = L 3 k = 1 -pl([r 2 k ]) 2 + pl([r 2 k ]) 3 , L 2 k = pl([r 2 k ]) 3 .
Having the same likelihood,

[x 1 k ] and [x 3 k ] are thus indistin- guishable. If furthermore pl([x 2 k ]
) is close to 1 then all the three likelihoods are quasi equal.

V. REINFORCED LIKELIHOOD BOX PARTICLE FILTER

A. Assumptions Consider system (Σ) with Assumptions (A), (B2), (D1) and (D2). Under these assumptions, Info (a) and Info (b) are concerned for likelihood computation.

Remark 4: The above measurement assumptions concern sensor errors (B2) and model (stochastic) uncertainties (D1). By (B2), the measurements are intervals

[y k ]. Regarding Remark 2, it is necessary to replace [y k ] ← [y k ] -[v k ] where [v k ] is the 99.7% confidence interval of v k determined by [v k ] = µ k -3 diag(Σ k ), µ k + 3 diag(Σ k ) , (8)
as proposed in [START_REF] Tran | Evidential box particle filter using belief function theory[END_REF]. This treatment generalizes Remark 2.

B. Method and Algorithm

The proposed method is named by Reinforced Likelihood Box Particle Filter (RLBPF) aiming to benefit advantages of existing criteria and LCMs and also attaining a gain in computation time. The method is based on the following two Propositions.

Proposition 1: Let {s i } i=1:n be a set of scores so that 1 ≤ s 1 ≤ s 2 ≤ ... ≤ s n and having an increasing order of preferability. Let Φ, Φ be normalizing functions so that (L 1 , ..., L n ) = Φ(s 1 , ..., s n ) = (s 1 , ..., s n )/ n k=1 s k and ( L1 , ..., Ln ) = Φ(s 1 , ..., s n ) = (e s1 , ..., e s2 )/ n k=1 e s k . Then the total distinguishability of the set { Li } i=1:n is greater than the total distinguishability of the set {L i } i=1:n .

Proof 1: One must prove i<j (s j -s i ) n k=1 s k < i<j (e sj -e si ) n k=1 e s k .
One can expresses:

LHS = n k=2 s k -(n -1)s 1 + n k=3 s k -(n -2)s 2 + n k=1 s k + ... + s n -s n-1 n k=1 s k = (n -1) -2 n-1 i=1 (n -i) s i n k=1
s k and the RHS has a similar expression form with replacing s i by e si , i = 1, ..., n. So :

RHS -LHS = 2 n-1 i=1 (n -i) s i n k=1 s k - e si n k=1 e s k = 2 n-1 i=1 (n -i) k =i (e s k s i -e si s k ) n k=1 s k n k=1 e s k = 2 n-1 i=1 n k=i+1 (k -i) (e s k s i -e si s k ) n k=1 s k n k=1 e s k For k > i, put ∆ ik = s k -s i , then e s k s i -e si s k = e si s i (e ∆ ik -1) -e si ∆ ik > e si ∆ ik (s i -1) ≥ 0
since e ∆ ik -1 > ∆ ik and using the assumption s i ≥ 1, ∀i = 1, ..., n. The proposition is proved. Proposition 2: Let Φ be a normalizing function so that Φ(a 1 , ..., a m ) = (a 1 , ..., a m )/ m k=1 a k , ∀{a i } i=1:m . Let {s i } i=1:n be a set of scores with an increasing order of preferability so that s 1 ≤ ... ≤ s n and {s i } i=k:n its truncated set with k ≥ 1. Let (L 1 , ..., L n ) = Φ(s 1 , ..., s n ) and ( L1 , ..., Ln-k+1 ) = Φ(s k , ..., s n ). Then

L i ≤ Li-k+1 , k ≤ i ≤ n, (9)
and the distinguishability of the truncated set is so that

L j -L i ≤ Lj-k+1 -Li-k+1 , k ≤ i < j ≤ n. ( 10
)
Proof 2: The proposition is proved by applying directly related definitions.

Principle: The scores {J 1 k , ..., J M k } are associated to M considered particles as follows:

J i k = d i k,1 + d i k,2 V i k , i ∈ {1, ..., M }, (11)
where

• d i 1 = d H ([y k ], [ŷ i k ]) (Hausdorff distance), • d i 2 = mid([y k ]) -mid([ŷ i k ]) 2 (Euclidean norm), • V i k = Vol([ŷ i k ] \ [y k ]) Vol([ŷ i k ]) = 1 - Vol([ŷ i k ] ∩ [y k ]) Vol([ŷ i k ])
, where Vol(.) is defined by (4). Thereby, J i k 's measure the closeness between [ŷ i k ]'s and [y k ] via both a kind of maximum distance d i k,1 and a kind of concentric tendency measure d i k,2 . J i k 's also take into account the size of intersections [ŷ i k ] ∩ [y k ] via the volume proportions V i k 's. Consequently, J i k 's exploit at the same time Info (a) and Info (b) and meets both Criterion 1 and Criterion 2. The particles having small scores are preferable, or equivalently the set {J i k } i=1:M has a decreasing order of preferability.

Then, the following strategy is applied: + {J i k } i=1:M is sorted in an ascending direction and a reduction percentage R% is applied, i.e. N h = (100 -R)%M particles corresponding to the scores {J i k } i=1:N h of the sorted set {J i k } i=1:M are retained. + the likelihoods (or weights) are computed as:

W i k = exp{-J i k + c} N h p=1 exp{-J i k + c} , i = 1 : N h , (12)
where c ≥ J N h k + 1 to let -J i k + c ≥ 1, i = 1 : N h . This strategy ensures an increasing of the distinguishability of the computed likelihoods thanks to Proposition 1 and 2.

After computing estimate [x k+1 ] according to (2), a backward estimate is added as follows:

[x k ] = hull{[x i k ]}, (13)
for those [x i k ]'s correspond to W i k 's just computed to reduce the conservatism due to interval computation.

Algorithm 1 Reinforced Likelihood Box Particle Filter

1: Initialization: [x 0 ] ≡ [x 0 ], R%, M , [u k ], [w k ], [y k ], [µ k ], [Σ k ], k = 1, ..., N , N h = (100 -R)%M . 2: for k = 1, 2, 3, ...N do 3: Partition [x k-1 ] to M disjoint boxes {[x i k-1 ]} i=1:M 4: [x i k ] = [f ]([x i k-1 ], [u k ], [w k ]) , i = 1, ..., M 5
:

[ŷ i k ] = [h]([x i k ], [u k ]) , i = 1, ..., M 6:
Compute J k = (J 1 k , ..., J M k ) using [START_REF] Nassreddine | State estimation using interval analysis and belief-function theory: Application to dynamic vehicle localization[END_REF].

7:

Sort J k in ascending direction:

J k = sort(J k ) 8:
Hold N h first values:

J k = J k (1 : N h ) 9: Compute W i k , i = 1 : N h , using (12)
.

10:

[x k ] = N h i=1 W i k .[x i k ]; [x k-1 ] = hull({[x i k ]} i=1:N h ) 11: end for
Remark 5: BPFs often use a non large number of particles to gain computation time and reduce the loss of a guaranteed estimation. Consequently, the resampling/repartition step happens at every or only after a few iterations. Thus, the fact that we hold previous weights and update them has no significant effect while this effect might not be quantified easily. Furthermore, conditions under which this procedure is implemented base usually on some heuristic choice of a threshold. It is also an issue of discussion but out of the scope of the present paper. Therefore, the RLBPF uses a reasonable small number of particles, performs the repartition at each iteration and strengthens the likelihood computation and the estimation by more efficient strategies.

VI. SIMULATION

Consider the following nonlinear system which was used as an illustration example in [START_REF] Tran | Evidential box particle filter using belief function theory[END_REF]. It is used in the present work to compare the proposed method (RLBPF) to the one (EBPF) in the reference, for the system states estimation.

x k+1 = α k,1 1 1 -α k,1 α k,2 x k + diag{β k,1 , β k,2 }u k + diag{20, 10}w k , y k = x k,2 x k /10 + v k , (14)
with

x k = (x k,1 , x k,2 ) T , u k ∈ [u k ] = 5([15, 17], [-7, -5]) T , w k ∈ [w k ] = ([-1, 1], [-1, 1
]) T /10 2 and for i ∈ {1, 2}:

β k,i = 0.5 + e T i x k /20 2 , α k,i = (0.2 + e T i δ k /20)(2β k,i -0.5), e 1 = (-1, 1) T , e 2 = (1, 2) T , δ k ∈ [δ k ] = 10[w k ]. The initial state is x 0 = [90, 80] T with [x 0 ] = ([85, 103], [75, 91]) T , the number of iteration N = 10 4 and v k ∼ N (µ k , Σ k ) where µ k ∈ [µ k ] = ([-1, 1], [-1, 1]) T and Σ k ∈ [Σ k ] = [90, 200].diag{1, 1}.
In order to evaluate how a LCM of a filter (e.g. RLBPF or EBPF) sharing Scheme 1 improves the efficiency of the estimation, we propose to compare the result of the filter with those of the following basic scenarios of Scheme 1:

• Scenario 1: Use the contraction step without partition (1 box); • Scenario 2: Use equi-likelihood 1 M without contraction step (M boxes); • Scenario 3: Use equi-likelihood 1 M with contraction step (M boxes). These basic scenarios can be considered as simple BPFs. The reason of this proposition is that, in some applications, using solely the contraction step, the algorithm performance has been rather good and the efficiency brought by the LCM might be insignificant. The same manner might happen for the other scenarios. In addition, following indexes, proposed in [START_REF] Tran | Evidential box particle filter using belief function theory[END_REF], are used for performance evaluations: where RM SE is the root mean squared error upper bound.

Let's consider the three basic scenarios previously defined. Table I shows that using only the contraction step gives no good performance results in terms of RM SE and E indexes (Scenario 1). Comparing Scenarios 2 and 3, it is shown that the contraction step brings a poor efficiency to the use of equilikelihood. Next, the simulation of RLBPF and EBPF applied to system (14) is considered. Since, in [START_REF] Tran | Evidential box particle filter using belief function theory[END_REF], the point value measurements y k are given, in this simulation we also use such an assumption and get [y k ] = y k -[v k ] where [v k ] is given by [START_REF] Freitas | Autonomous crowds tracking with box particle filtering and convolution particle filtering[END_REF]. The reduction percentage R = 20% is applied for RLBPF and the particle number M = 9 is applied for both methods. Table II shows the better performance of RLBPF versus EBPF in terms of RM SE, E indexes and the computation time (reduced more than 60%). In addition, EBPF performance is not better than those of Scenarios 2 and 3 in all indexes and computation time. In contrast, the RM SE and E indexes provided by RLBPF are better than those of Scenario 2 and 3, while its computation time is a compromise between those of the two Scenarios. A general scheme is provided to generalize the specificity of BPFs. The likelihood computation methodology is investigated. This analysis point out the disadvantages of existing filters and opens a way to improve the computed likelihoods by making them more reliable using a reinforcement method. A strategy is proposed to evaluate the performance of this class of filters. The simulation highlights the efficiency of the RLBPF in gain of computation time and evaluation indexes.

Fig. 1 .

 1 Fig. 1. Likelihood computation methodology diagram

2 ,

 2 RM SEj = sup N k=1 (x k,j -[x k,j ]) 2 /N 1/Ej = N k=1 width([x k,j ])/N , Oj = N k=1 1(x k,j ∈ [x k,j ])/N 100%, j = 1, ..., nx,

Fig. 2 .

 2 Fig. 2. RLBPF versus EBPF

  B3) The measurements are assumed to be accurate in the sense that [y k ] h(x k , u k ) ≡ h(x k , u k , 0) (the zero noise case), where x k is the real state. Assumption (C): Measurement Stochastic Uncertainty (C1) v k are additive noises with known density p v . (C2) The measurements are point values y k . Assumption (D): Measurement Mixed Uncertainty (D1) v k are additive Gaussian noises with unknown mean µ k ∈ R ny and covariance

nw state dynamic disturbance and v k ∈ R ny measurement noise. Assumption (A): State Process Uncertainty u k and w k are unknown and belong to known intervals [u k ] and [w k ] respectively. Assumption (B): Measurement Bounded Uncertainty (B1) v k (unknown) belongs to known interval [v k ]. (B2) The measurements are intervals [y k ].

(