Vahid Khanagha
email: vahid.khanagha@inria.fr

Khalid Daoudi
email: khalid.daoudi@inria.fr

Oriol Pont
email: oriol.pont@inria.fr

Hussein Yahia
email: hussein.yahia@inria.fr

Phonetic segmentation of speech signal using local singularity analysis

Keywords: Nonlinear speech processing, Multiscale signal processing, phonetic segmentation of speech signals

This paper presents the application of a radically novel approach, called the Microcanonical Multiscale Formalism (MMF) to speech analysis. MMF is based on precise estimation of local scaling parameters that describe the inter-scale correlations at each point in the signal domain and provides efficient means for studying local non-linear dynamics of complex signals. In this paper we introduce an efficient way for estimation of these parameters and then, we show that they convey relevant information about local dynamics of the speech signal that can be used for the task of phonetic segmentation. We thus develop a two-stage segmentation algorithm: for the first step, we introduce a new dynamic programming technique to efficiently generate an initial list of phoneme-boundary candidates and in the second step, we use hypothesis testing to refine the initial list of candidates. We present extensive experiments on the full TIMIT database. The results show that our algorithm is significantly more accurate than state-of-the-art ones.

Introduction

The existence of highly non-linear and turbulent phenomena in production process of the speech signal has been theoretically and experimentally established [START_REF] Kaiser | Some observations on vocal tract operation from a fluid flow point of view[END_REF][START_REF] Teager | Evidence for nonlinear sound production mechanisms in the vocal tract[END_REF]. However, most of the classical approaches in speech processing are based on linear techniques that may not adequately capture the complex dynamics of the speech signal. Some of the well-known examples of such non-linear phenomena include: the existence of turbulent sound source in production process of unvoiced sounds [START_REF] Kaiser | Some observations on vocal tract operation from a fluid flow point of view[END_REF][START_REF] Kubin | Speech coding and synthesis[END_REF][START_REF] Maragos | Fractal dimensions of speech sounds: Computation and application to automatic speech recognition[END_REF], the existence of a time spread and turbulent component for the excitation source of plosives [START_REF] Quatieri | Discret-time speech signal processing: principles and practice[END_REF] (that is idealized as an impulse in the linear framework) and the evidences regarding characterization of voiced sounds by highly complex air flows like jets and vortices [START_REF] Teager | Evidence for nonlinear sound production mechanisms in the vocal tract[END_REF].

There is a broad spectrum of non-linear speech analysis approaches [START_REF] Kubin | Speech coding and synthesis[END_REF][START_REF] Fandez-Zanuy | Nonlinear speech processing: Overview and applications[END_REF][START_REF] Little | Biomechanically informed nonlinear speech signal processing[END_REF][START_REF] Mclaughlin | Nonlinear methods for speech analysis and synthesis[END_REF]. One common approach is to consider speech as a realization of a non-linear dynamical system and use the tools and methods in the study of such systems to catch the non-linear character of the speech signal. In these approaches, the non-linear aspects of the speech signal are associated to the turbulent nature of airflow in the vocal tract [START_REF] Kaiser | Some observations on vocal tract operation from a fluid flow point of view[END_REF][START_REF] Maragos | Fractal dimensions of speech sounds: Computation and application to automatic speech recognition[END_REF][START_REF] Teager | Evidence for nonlinear sound production mechanisms in the vocal tract[END_REF][START_REF] Little | Exploiting nonlinear recurrence and fractal scaling properties for voice disorder detection[END_REF], which in turn justifies the use of methods from the study of turbulent systems for speech analysis [START_REF] Fandez-Zanuy | Nonlinear speech processing: Overview and applications[END_REF].

For instance in [START_REF] Pitsikalis | Analysis and classification of speech signals by generalized fractal dimension features[END_REF], it is assumed that the finite-dimensional embedded attractor of speech in the phase-space shares common features with its unknown non-linear production system. Consequently, topological features of the unknown dynamical system can be estimated from its embedded attractor. Correlation dimension and general dimension are two examples of such invariant quantities that are employed in [START_REF] Pitsikalis | Analysis and classification of speech signals by generalized fractal dimension features[END_REF] to define new features for speech analysis. Another characteristic of a dynamical system which remains intact by embedding procedure is Lyapanov Exponents (LE). LEs characterize the global degree of chaos or predictability in a dynamical system. They are shown in [START_REF] Kokkinos | Nonlinear speech analysis using models for chaotic systems[END_REF] to be useful features for phoneme classification. Also, a positive value of LE (which can be considered as an indicator for chaoticity) is observed for most of the phonemes.

These experiments show the potential of employing the tools and concepts coming from the study of chaotic, turbulent systems for speech analysis. However, most of these approaches suffer from practical and theoretical limitations. Apart from the rather involved mathematical procedures, the major limitation of most of these methods is their underlying assumption of stationarity, which is clearly not valid for a rapidly time-varying signal like speech. Consequently, these tools and methods can only be applied on isolated parts of the signal that are assumed stationary, like isolated phonemes. However, another issue appears when dealing with isolated phonemes: the length of a single phoneme is usually less than sufficient for precise estimation of parameters like LEs. This is why sustained vowels are used in [START_REF] Banbrook | Speech characterization and synthesis by nonlinear methods[END_REF] instead of naturally uttered ones. In another study [START_REF] Banbrook | Speech characterization and synthesis by nonlinear methods[END_REF] extensive experiments are conducted so as to choose the optimal dynamical model that can be used for relatively more precise estimation of LEs from very short signals like phonemes (to avoid the use of sustained vowels). Other than such practical issues, the fact that these approaches are based on the estimation of some global parameters, limits their usability for accessing many of the important local dynamical properties of the speech signal. These methods can generally recognize the existence of highly complex structures, but they can not provide local access to these structures.

In this paper, we use a radically novel framework coming from the study of highly disordered turbulent flows in statistical physics that does not suffer from the aforementioned limitations. This formalism offers new ways to evaluate quantitatively and locally the degree of predictability at every point in the signal domain. This is done trough the accurate estimation of local scaling parameters called Singularity Exponents (SE), which not only allows the identification of complexity inside signals but also allows to localize geometrically, inside the signal, where the complexity happens and how it organizes itself. In this paper we present a novel efficient method in estimation of SEs which is directly related to the degree of predictability around a point. We show how they convey meaningful information about some local dynamics of speech. Then, to demonstrate the usefulness of these exponents, we show how they can be easily used to develop an accurate and efficient phonetic segmentation algorithm.

Given that most of the recent developments in speech technology strongly rely on corpusbased methodologies that require the availability of precisely time-aligned labels of phonemes, phonetic segmentation is a very important task. The most precise segmentation method is to manually perform the task, which is extremely complex and expensive. Hence, automatic phonetic segmentation is of great importance and interest. The most frequent approach is to use an HMMbased phonetic recognizer supplied with the associated phonetic transcription and performing a forced alignment [START_REF] Torre-Toledano | Automatic phonetic segmentation[END_REF] to detect phoneme boundaries. Although such methods provide very good segmentation performances, they suffer from imposing linguistic constraints to the segmentation algorithm. This makes the algorithms restricted to the database used for training and hence, they can not be applied to different databases of different languages, contexts or accents. Moreover, it is not always possible to access phonetic transcriptions to perform forced time-alignment. Such HMM based methods, or any other approach which relies on an externally supplied transcription of the sentence for determining phoneme boundaries are called Text-Dependent (TD) segmentation methods. Text Independent (TI) methods on the other hand rely exclusively on the acoustic information contained in the speech signal itself and do not require externally supplied transcriptions. TI methods can be used to simplify the manual labeling process for development of speech corpora. They may play a very important role in phoneme-based multi-lingual speech recognition systems [START_REF] Esposito | Text independent methods for speech segmentation[END_REF].

The method we use in this paper belongs to the TI class of segmentation methods. We try to exploit our precisely estimated SEs to identify the phoneme boundaries. We employ a classical two-step methodology in change detection to identify the points in which local properties of SEs are changing. First, we use the measure we introduced in [START_REF] Khanagha | A novel text-independent phonetic segmentation algorithm based on the microcanonical multiscale formalism[END_REF], which represents the instantaneous average of SEs to select a set of candidates for phonetic boundaries. To do so, we use the piecewise-linear approximation of this measure applied to signal itself and its low-passed filtered version. In the second step, the initial list of candidates is refined using dynamic windowing and Log-Likelihood Ratio Test (LLRT). This two-step approach is similar to traditional segmentation methods that are based on a boundary pre-selection step followed by statistical tests to make the final decision [START_REF] Andre-Obrecht | A new statistical approach for the automatic segmentation of continuous speech signals[END_REF][START_REF] Almpanidis | Robust detection of phone boundaries using model selection criteria with few observations[END_REF]. Some of the preliminary results of our algorithm are presented in [START_REF] Khanagha | Improving text-independent phonetic segmentation based on the microcanonical multiscale formalism[END_REF]. In this paper, we describe the detailed procedure for SE estimation (a novel procedure that we have particularly adopted for phonetic segmentation problem (section 2.1)). Also, our novel dynamic-programming-based technique for piece-wise-linear approximation [which has crucial role in efficiency and text-independency of our segmentation algorithm] is fully presented in section 4.2. Moreover, we present more extensive experimental results (the test part of TIMIT database is included).

We evaluate our method over full TIMIT database including its train and test sets and we report all of the well-known performance measures such as detection, insertion and over-segmentation rates along with two global measures which simultaneously take these three into account. By doing so for several degrees of error tolerance, we provide results which are easy to interpret and to compare with. We compare our results with two state-of-the-art methods [START_REF] Dusan | On the relation between maximum spectral transition positions and phone boundaries[END_REF][START_REF] Esposito | Text independent methods for speech segmentation[END_REF] which have reported their results on these known datasets (train and test dataset respectively). [START_REF] Qiao | Unsupervised optimal phoneme segmentation: Objectives, algorithm and comparisons[END_REF] and [START_REF] Lin | Phone boundary detection using sample-based acoustic parameters[END_REF] also report on the train dataset but the algorithms they propose are not fully unsupervised. The former assumes prior knowledge of the number of phonemes in the utterance, while the latter uses prior knowledge of all manual transcriptions to train a neural network (and hence the "textindependence" property of the system is compromised). Applied to the train set, our algorithm significantly outperforms [START_REF] Dusan | On the relation between maximum spectral transition positions and phone boundaries[END_REF]. Consequently, the application of the algorithm to the test dataset, without any tuning of parameters, shows comparable results w.r.t the best performance of [START_REF] Esposito | Text independent methods for speech segmentation[END_REF] which is achieved after extensive tuning of parameters on the whole dataset.

The paper is structured as follows: in Section 2 we briefly introduce the basic concepts of the MMF and the detailed procedures for computation of SEs. In section 3 some examples about the informativeness of SEs are presented and the detailed algorithm for using them in the TI phonetic segmentation task is presented in section 4. The experimental results are presented in Section 5 and we draw our conclusion in Section 6.

The Microcanonical Multiscale Formalism

The approach we use in this paper for nonlinear speech analysis is based on the Microcanonical Multiscale Formalism (MMF) [START_REF] Turiel | Microcanonical multifractal formalism: a geometrical approach to multifractal systems. part 1: singularity analysis[END_REF], which allows the study of local geometrico-statistical properties of complex signals from a multi-scale perspective. The particularity of this formalism compared to its older counterparts (the canonical formulations [START_REF] Frisch | Turbulence: The legacy of A.N. Kolmogorov[END_REF]) is that it considers local quantities defined at each point of the signal domain, instead of global quantities like moments and structure functions. Central to the formalism is the computation of Singularity Exponents (SE) at each point in the signal domain. When correctly defined and estimated, these exponents alone can provide valuable information about local dynamics of complex signals and have been successfully used in many applications ranging from signal compression to inference and prediction [START_REF] Grazzini | Edge-preserving smoothing of high-resolution images with a partial multifractal reconstruction scheme[END_REF][START_REF] Yahia | Motion analysis in oceanographic satellite images using multiscale methods and the energy cascade[END_REF][START_REF] Maji | Inferring information across scales in acquired complex signals[END_REF][START_REF] Maji | Towards multiscale reconstruction of perturbated phase from hartmann-shack acquisitions[END_REF][START_REF] Badri | Computer graphics effects from the framework of reconstructible systems[END_REF]. SE values can be used to analyze the local scaling behavior of complex signals. It is shown [START_REF] Turiel | Microcanonical multifractal formalism: a geometrical approach to multifractal systems. part 1: singularity analysis[END_REF] that for a given time-instant t, the smaller the value of SE is, the higher unpredictability is at this point. This property can be used to identify an important subset of signal samples that carry most of the information content of the signal: the Most Singular Manifold (MSM). MSM is defined as the subset of signal samples whose SE values are smaller than a threshold. It has been established that the critical transitions of the system occur at these points and this fact has been successfully used in many applications [START_REF] Turiel | Microcanonical multifractal formalism: a geometrical approach to multifractal systems. part 1: singularity analysis[END_REF][START_REF] Yahia | Motion analysis in oceanographic satellite images using multiscale methods and the energy cascade[END_REF].

The Singularity Exponent (SE), h(t), for a signal s(t), can be estimated by evaluation of the power-law scaling behavior of a multi-scale functional Γ r over a set of fine scales r:

Γ r (s(t)) = α(t) r h(t) + o r d+h(t) r → 0 (1)
where Γ r (•) can be any multi-scale functional complying with this power-law. The term o r h (t) means that for small scales the additive terms are negligible and thus h(t) dominantly quantifies the multi-scale behavior of the signal at time t. α(t) is a quantity that is independent of the scale and can be separated from h(t).

An important aspect in the MMF is the definition Γ r (•) such that the inter-scale power-law correlations in form of Eq. ( 1) are revealed. As suggested in [START_REF] Turiel | Microcanonical multifractal formalism: a geometrical approach to multifractal systems. part 1: singularity analysis[END_REF], one reasonable choice is to integrate all the variations (dµ D ) of the signal on a ball B r of radius r centered at t:

Γ µ r (s(t)) = B r dµ D (2)
while the variations being defined as:

dµ D = Ds dr ( 3
)
where D is an appropriate differential operator such as the norm of the gradient ∇s , whose choice is motivated in [START_REF] Turiel | Microcanonical multifractal formalism: a geometrical approach to multifractal systems. part 1: singularity analysis[END_REF] by the special power spectrum scaling of natural images. The use of the norm of the gradient helps to avoid the domination of the integral in Eq. ( 2) by the mean value of s, which may be far from zero due to the lack of stationarity. The resulting multi-scale measure is related to the typical characterization of intermittency in turbulence and describes the kinetic energy dissipation at scale r on B r for a turbulent velocity field [START_REF] Frisch | Turbulence: The legacy of A.N. Kolmogorov[END_REF][START_REF] Turiel | The multi-fractal structure of contrast changes in natural images: from sharp edges to textures[END_REF].

Estimation of the singularity exponents

Ideally, if the power-law scaling behavior of Eq. ( 1) holds for a given measure Γ r (•), the measurements at different scales are directly related. In this case, SEs can be estimated by evaluation of Eq. (1) using only one measurement at the finest scale: if one takes the logarithm of both sides of Eq. (1), α(t) appears as an additive term that would be negligible for an adequately small scale r 0 [START_REF] Turiel | Numerical methods for the estimation of multifractal singularity spectra on sampled data: A comparative study[END_REF]. Now if we drop the constant shift d, one can have an estimate of the SE as:

h(t) = log Γ µr 0 (s(t)) Γ µr 0 (s(•)) log r 0 ( 4
)
where Γ µ r 0 (s(•)) is the average of the functional over the whole domain of the signal and serves to diminish the additive factor due to α(t), as explained in [START_REF] Turiel | Numerical methods for the estimation of multifractal singularity spectra on sampled data: A comparative study[END_REF]. There, it is shown that the additive factor due to α(t) is negligible provided that the finest accessible scale (r 0 ) is small enough. By using the punctual estimation one preserves the finest accessible resolution in SE estimation. However, in practice we do not usually have access to the physical finest scale (due to discretization limitations) and hence the resulting estimate (Γ µ r 0 ) might be unstable and perturbed by the additive term α(t). This perturbation can be corrected by incorporating the measurement at coarser scales. But this has the risk of compromising the resolution and the measurements at coarser sampling scales may simply combine the information carried by adjacent samples. For the application we consider in this paper, resolution is of very high importance. So we only add the measurement at the second finest scale (Γ µ r 1 ) to regularize Γ µ r 0 as:

Γ κ r (s(t)) = κ t Γ µ r 0 (s(t)) (5)
where the regularization term κ t is defined as the quotient of the measurements at two scales:

κ t = Γ µ r 1 (s(t)) Γ µ r 0 (s(t)) (6)
where r 1 is the next coarser scale (r 1 > r 0 ). In the ideal case both measurements are completely correlated and κ t = 1. But in practice, κ t serves to compensate for the perturbations of a possibly noisy measurement made at the finest scale r 0 . We thus replace Γ µ r 0 in Eq. ( 4) with Γ κ r (s(t)) to compute h(t). A summary of the overall procedure for computing the SE is provided in Alg. 1. Note that the procedure is written for a 1-D discrete time signal s[n] sampled from s(t) with the sampling frequency of f s. For such a discrete time signal, the finest accessible scale is the one provided by the sampling frequency. Consequently an estimate of the sum of variations at the finest scale (Eq. ( 2)) would simply be Γ µ r 0 (s[n]) = |s[n]s[n -1]| and at the second finest scale, we have two measurements to sum up as

Γ µ r 1 (s[n]) = |s[n] -s[n -1]| + |s[n + 1] -s[n]|.
Algorithm 1 The procedure for computation of the SEs for a discrete time signal

1: Γ µ r 0 (s[n]) ← |s[n] -s[n -1]| 2: Γ µ r 1 (s[n]) ← |s[n] -s[n -1]| + |s[n + 1] -s[n]| 3: κ t [n] = Γ µr 1 (s[n]) Γ µr 0 (s[n]) 4: Γ κ r (s[n]) = κ t [n]Γ µ r 0 (s[n]) 5: h[n] = log Γ κ r (s[n]) Γµ r 0 (s[n]) log r 0 6: ⊲ r 0 = 1 f s ,
• denotes the average value over the whole sentence

A singularity analysis of phonemes

In this section we present some preliminary observations on how the temporal evolution of h[n] conveys instructive information that can be used to detect phoneme boundaries. In Fig. 1-top a speech signal from the TIMIT database is shown along with the reference phoneme boundaries extracted from the manual transcription of TIMIT. The inter-phoneme variability of exponents can be further demonstrated using the concept of MSM. As mentioned in section 2, MSM is defined as the set of samples in the signal domain whose value of h[n] is smaller than a threshold and it theoretically represents the most informative subset of samples in the signal domain. In our experiment, we set a global threshold h th for the whole speech signal and in each window of length 32ms around t, we count the number of samples with h[n] < h th . The result is shown in Fig. 1-bottom. It can be seen that the counter value is almost constant inside each phoneme (if we neglect the transition on phoneme boundaries where the sliding window partially covers both of the neighboring phonemes) and this constant value is different for adjacent phonemes. These observations suggest that SE do indeed convey valuable information about local dynamics of the speech, and in this very simple form can be readily used for phoneme segmentation. We next develop new tools to employ these properties for practical development of an automatic phonetic segmentation algorithm.

Application of singularity analysis to phonetic speech segmentation

We use a two steps procedure to employ the properties of SEs mentioned in section 3 for the task of automatic phonetic segmentation. We first define a measure that captures the interphoneme changes in the properties of SEs in a very simple form and use it to make a preliminary list of phoneme boundary candidates. Then we perform a dynamical windowing with a hypothesis testing to decide whether each candidate truly represents a phoneme transition.

The Accumulative function ACC

The simplest parameter that can quantify the changes in distribution of SEs between neighboring phonemes is their average. We consider h[n] as a random variable whose average is changing between adjacent phonemes. We thus search for the locations of changes in local averages of h[n] and consider them as the candidate phoneme boundaries. Local averages of this random variable can be estimated in a small sliding window. However, as SE estimations are available at the finest resolution and we are interested in preserving this resolution for the task of phonetic segmentation, we try to avoid windowing for estimation of averages. Instead, we use the running mean of h[n] defined as:

ACC[n] = n k=1 (h[k] -h) ( 7
)
where h is the global average of exponents (over the whole sentence) and N is the total number of speech samples. Indeed, inside the boundaries of each phoneme, the slope of ACC[n] would be an estimate for the local average of the h[n] minus the constant term h. The resulting functional is plotted in Fig. 2. Just as we expected, this new functional reveals the changes in distribution in a more precise way. Indeed ACC[n] is almost linear inside each phoneme, while there is a clear change in the slope at the phoneme boundaries. Extensive observations over different sentences confirm this behavior and show that these slope changes are even happening at the boundaries between extremely short phonemes, such as stops. Consequently, detection of the locations of changes in the slope of the ACC can lead to an automatic phonetic segmentation algorithm. A very simple solution is to fit a piece-wise-linear curve to the ACC and take the break-points as the candidates of change in distribution of SE, i.e. phoneme boundary candidates.

Piece-wise-linear approximation of ACC

For a speech signal of length N that contains K phoneme boundaries, we aim at fitting a piece-wise-linear curve to the corresponding ACC [n] such that the breakpoints are located on the phoneme boundaries that we are searching for. If we denote the m-piece linear approximation of ACC[n] (n ∈ [1, N]) by LA m (ACC, 1, N) and show the corresponding mean squared error of this approximation by E 1→N m the following optimization problem can be used to reach our goal (conventional method):

• find LA m (ACC, 1, N) such that m is minimized and E 1→N m < ǫ.
where ǫ is a reasonable threshold that is conventionally set by trial and error. Ideally, we wish that m = K so that there are no false detections or deletion of boundaries. An intuitive solution for this problem is to recursively search through all the points in the domain of ACC[n] to find the breakpoints (one after another) such that at each step the least possible value of approximation error is achieved. This method is addressed in [START_REF] Little | Generalized methods and solvers for noise removal from piecewise constant signals: Parts i and ii[END_REF] for denoising of piece-wise constant functions and it is argued that when K << N, such a greedy search may result in a more efficient and more accurate approximation compared to other solvers. However, considering the large values of N and K in the particular problem of phonetic segmentation, the computational complexity of this method (which increases with K) is still unbearable when dealing with large databases of the speech signal. Also, since we don't know the number of phoneme boundaries (K), the stop condition of the algorithm (the value of ǫ) is vague. We thus develop a more efficient method whose computational load is of K.

Our optimization method is motivated by the dynamic programming approach of [START_REF] Ihm | Scientific visualization of physical phenomena[END_REF] where the above optimization problem is solved by recursive solution of easier sub-problems of the same type. In our method, instead of recursive search for K breakpoints over the whole signal each time, we perform a left-to-right search that passes through the signal only once: from left to right, each time a break-point is selected, the part of ACC before that breakpoint is removed from the search in the next iterations. More specifically, we start from n = 1 and search for the first point n 1 such that the error of the 1-piece approximation of ACC[n] with a straight line from ACC [START_REF] Kaiser | Some observations on vocal tract operation from a fluid flow point of view[END_REF] to ACC[c i ] becomes larger than a threshold ǫ. Clearly the location of n 1 is dependent on the choice of ǫ, which is not desirable because we are searching for specific fixed location of phoneme boundaries. We overcome this issue in the second step: we first opt for a relatively large threshold to ensure that noise-like variations are not mistaken as breakpoints. Once n 1 is found, we search for the breakpoint c 1 ∈ [1, n 1 ] that once taken as a breakpoint, the 2-piece linear approximation 1 error in this segment is minimized. As such, the breakpoint c 1 is found through a local minimization procedure and not just a thresholding procedure and hence, we expect the overall algorithm to be less sensitive to the choice of the threshold ǫ. Once c 1 is fixed, we perform the exact same procedure, this time starting from c 1 to find the second breakpoint and repeat the procedure until the last sample is reached. This procedure is detailed in Alg. 2. ¢ Algorithm 2 PLA procedure

1: i ← 1; c 0 ← 1; k 1 ← 2; 2: while k 1 < N do 3: if E c i-1 →k 1 1 > ǫ then; 4: n i ← k 1 ; 5: for k 2 ∈ (c i-1 • • • n i ) do
6:

E c i-1 →n i 2 [k 2 ] ← E c i-1 +1→k 2 1 + E k 2 +1→k 1 1 ; 7:
end for 8:

i ← i + 1;

9:

c i = argmin k 2 E c i-1 →n i 2 [k 2 ]; 10: k 1 ← c i ; 11:
end if ⊲ At each iteration c i-1 denotes the previously taken phoneme candidate.

We emphasize that at each iteration, c i is determined through a greedy minimization (Alg. 2 line 9) and not through a thresholding operation. This greedy minimization decreases the impact of the threshold ǫ (Alg. 2 line 3). We will later show in the experimental results that the algorithm is not so sensitive to the selection of ǫ. The insensitivity to the selection of algorithm parameters is an important property for text-independent phonetic segmentation, where no data is available to train the parameters.

The two-step segmentation algorithm

The piece-wise linear reconstruction of ACC as explained in section 4.2 is readily a working phonetic segmentation algorithm. However, by performing detailed analysis on the performance of this algorithm we found that it can not fully exploit the strength of SE for phone boundary detection. Motivated by some observations on the behavior of this simple algorithm at some particular phoneme transitions, we propose a two-step algorithm where we first pre-select candidate boundaries and then use statistical hypothesis test to make the final decision.

We first observed that some of the missed boundaries by the method of section 4.2 correspond to transitions between fricatives/stops and vowels. We also observed that transitions between speech and low energy segments (such as pauses and epenthetic silence) display strong and easy to detect changes in the slopes of ACC. Indeed, as shown in Fig. 1-middle, the SEs of low energy segments have high positive values, while they are mostly negative in active speech segments. Motivated by these observations and the fact that fricatives/stops are essentially high-band signals, we decided to compute ACC on a low-pass filtered version of the utterance. By doing so fricatives/stops-vowels transitions will be converted into silence-speech transitions which are easier to detect as shown in Fig. 3. It is known that most of the spectral energy of fricatives is located above 2000Hz and, for most stops, the active frequency bands start at 1800Hz [START_REF] Ladefoged | A Course in Phonetics, Chapter 8: Acoustic[END_REF]. We thus choose the cut-off frequency of the low-pass filter as 1800Hz. The second and most important observation we made is related to the statistical distribution of SEs. We observed that some missed boundaries correspond to neighboring phonemes which although have a quite strong difference in their SE distribution, the change in their averages is not strong enough to be reflected as a change in the slope of ACC and thus it is not captured by the simple linear approximation procedure. It is then natural to consider a statistical hypothesis test over SE distributions in our segmentation algorithm in order to detect such boundaries.

Motivated by these observations, we use a similar approach as in [START_REF] Almpanidis | Robust detection of phone boundaries using model selection criteria with few observations[END_REF] develop a new segmentation algorithm which consists of 2 steps. First, we apply Alg. 2 to the original signal and to its filtered version. We gather all the detected breakpoints and consider them as candidate boundaries. In the second step, we make the final decision on each candidate by performing a dynamic windowing over them followed by Log Likelihood Ratio Test (LLRT) over SE distributions of the original signal. We use a Gaussian hypothesis because our goal is to detect changes in mean and variance of SEs. More precisely, for each candidate c i we consider the large window Z = [c i-1 , c i+1 ] and the two smaller windows X = [c i-1 , c i ] and Y = [c i , c i+1 ]. We then compute LLR statistic to decide between the two hypotheses:

• H 0 : SE of Z are generated by a single Gaussian.

• H 1 : SE of Z are generated by two Gaussians on X and Y.

If H 1 is significantly likelier than H 0 , we select c i as a boundary. Otherwise, c i is removed from the candidates list. We emphasize here that SE of the filtered signal are used only in the first step. The final decision is exclusively made upon the information conveyed by SE of the original signal. We also emphasize that this new algorithm is still simple and efficient as the original one.

Experimental results

Our evaluation is carried out on the full Train and Test sets of the TIMIT database which respectively contain 4620 and 1200 sentences uttered by 462 and 120 speakers. We first develop our algorithm on a small dataset of 30 sentences extracted from the Train set. First we test the performance on the whole Train set and we compare the results with the state-of-the art methods. Then, to compare with [START_REF] Esposito | Text independent methods for speech segmentation[END_REF], and to show that the promising result on Train dataset is not the result of over-trained parameters, we evaluate it on the Test dataset which is balanced for the phonological coverage.

Performance measures

The segmentation quality can be evaluated and analyzed using three "partial" scores: the Hit Rate (HR), the False Alarm Rate (FA) and the Over Segmentation Rate (OS ). These three scores are defined in Table 1.

Table 1: Partial performance measures. N T is the total number of detected boundaries (correct and incorrect), N H is the number of correctly detected boundaries and N R is the total number of boundaries in the reference transcription

Measure

Formula Description

Hit Rate (HR)

HR = N H N R
The percentage of reference boundaries that are correctly detected Segmentation Rate (OS )

OS = N T -N R N R
OS shows how much more (or less) is total number of algorithm detections, compared to the total number of reference boundaries taken from the manual transcription.

False Alarm Rate (FA)

FA = N T -N H N T
The percentage of algorithm detections that are incorrect.

In order to assess the overall quality of a segmentation method, a global measure which simultaneously takes these scores in to account is required. A well known measure is the F 1value:

F 1 = 2 × PCR × HR PCR + HR ( 8
)
where PCR = 1 -FA. Another global measure, called the R-value has been recently proposed in [START_REF] Rasanen | An improved speech segmentation quality measure: the R-value, Proceedings of INTERSPEECH[END_REF]. This measure makes more emphasize on over-segmentation by arguing that better hit rates might be achieved by simply adding random boundaries without any algorithmic improvement. This measure evaluates how close one is to the ideal segmentation R = 1:

R = 1- |r 1 | + |r 2 | 2 r 1 = (1 -HR) 2 + OS 2 , r 2 = HR -OS -1 √ 2 (9)

Results : Train dataset

Using different sizes of tolerance windows, we provide comparison of segmentation results for 3 methods. In the first one, we give the results reported in [START_REF] Dusan | On the relation between maximum spectral transition positions and phone boundaries[END_REF]. We mention here that [START_REF] Dusan | On the relation between maximum spectral transition positions and phone boundaries[END_REF] report scores with 0ms, 10ms and 20ms tolerance windows. However, their approach is framebased with a 10ms frame step size and they convert each manual boundary to the closest frame position. Thus, 5ms has to be added to their window size in order to make a fair comparison with our sample-based approach, which has the finest possible resolution. In the second one, we provide the results using our original SE-based algorithm summarized in section 4.2, we call it SE-ACC. In the third one, we present the results obtained using our new algorithm described in section 4.3, we call it SE-LLRT. Table 2 presents HR, FA and OS for the 3 methods. The first observation is that SE-LLRT outperforms SE-ACC for the 3 scores and all tolerance windows. In particular, a significant improvement is made in FA and OS. This shows that, as expected, some of the insertions introduced by the curve fitting procedure have been corrected by the LLRT. The second observation is that SE-LLRT yields higher accuracy than [START_REF] Dusan | On the relation between maximum spectral transition positions and phone boundaries[END_REF]. In particular, the smaller the tolerance window, the higher the relative improvement is. This shows that SE-LLRT is better suited for high precision detection of phoneme boundaries. To this regard, we can mention another interesting comparison with [START_REF] Lin | Phone boundary detection using sample-based acoustic parameters[END_REF] which is also a sample-based segmentation method as ours. In [START_REF] Lin | Phone boundary detection using sample-based acoustic parameters[END_REF], it is reported that 43.5% of their 86.8% detection output is located within the first bin of the cumulative histogram of distances from true boundaries. This corresponds to 43.5% × 86.8%=37.75% hit rate with 7.5ms tolerance. With SE-LLRT we obtain 44% hit rate which is significantly more accurate. More importantly, the algorithm in [START_REF] Lin | Phone boundary detection using sample-based acoustic parameters[END_REF] is supervised (all manual transcriptions are used to train a neural network) while ours is fully unsupervised. 3 presents the performance of each of the 3 methods when evaluated using the global measures F 1 and R. The same observations we made above, still hold for the global performance evaluation. Indeed, SE-LLRT still outperforms SE-ACC for both F 1 and R. Moreover, about 6% (resp. 10%) improvement in R-value and 4% (resp. 10%) in F 1 -value is achieved for 25ms of tolerance (resp. 5ms and 15ms). This is a significant gain in accuracy that shows the strength of SEs in revealing the transitions fronts between phonemes.

Results : Test dataset

To the best of our knowledge, no detailed result is reported in the literature for the performance of a TI method on the TEST dataset. However, to facilitate later comparisons, we have provided the full segmentation results for this dataset, in the last columns of table 2 and 3. In [START_REF] Esposito | Text independent methods for speech segmentation[END_REF] the result for Test dataset is reported only for 20ms tolerance, as shown in table 4. It can be seen that the hit-rate of [START_REF] Esposito | Text independent methods for speech segmentation[END_REF] (8-Mel-bank coding scheme) is about 10% higher than that of ours, while its over-segmentation is about 13% worse. Hence a direct comparison in terms of hit-rate and over-segmentation is not meaningful. However, in terms of R-value the results are very close.

It should be mentioned that the results of [START_REF] Esposito | Text independent methods for speech segmentation[END_REF] are attained by performing an extensive parameter-tuning procedure over the whole Test dataset. Quite the contrary, in our case we have applied our algorithm to the Test dataset without any parameter tuning. Our intention was to ensure that the promising performance of the algorithm is not the result of over-training of parameters. Even more, we emphasize an important feature of our algorithm which is its insensitivity to the threshold of the linear curve fitting. Recalling that this threshold is applied to the mean-squared approximation error of the ACC, an initial value for the threshold is selected as 0.001 times the maximum value of ACC (which is normalized to unity). Fig. 4 shows how R-value varied for different thresholds. We used a subset of 30 randomly selected sentences to compute these values. It can be seen that with about 400% change of the threshold value, the variance of changes in R-value is less than 0.5%. Thus, it seems fair to claim that our algorithm is threshold-independent: SE-ACC does not have any other parameter to tune and for the SE-LLRT, the threshold for LLRT is the only parameter to tune. This is a major advantage as most of the TI methods require accurate threshold tuning. For instance, the results of [START_REF] Esposito | Text independent methods for speech segmentation[END_REF] are obtained by extensive optimization of three important parameters.

Conclusions

In this paper we used a novel framework in complex-signal analysis, called the MMF, to analyze the local dynamics of speech. We presented an efficient method for estimation of the singularity exponents which are the relevant features in this framework to identify critical transitions. Observing that temporal evolution of these exponents convey useful information about phoneme boundaries, we defined an accumulative functional upon these exponents which permits automatic detection of phone boundaries. This accumulative functional exhibits a piecewise-linear behavior with distinctive breakpoints at phoneme boundaries. We then developed a simple and efficient algorithm using dynamic programming to detect these break-points as the phoneme-boundary candidates (with special care for reducing computational complexity and also dependency on parameter tuning, which is necessary for text-independent phonetic segmentation). The final decision was made for each candidate, by performing a hypothesis test on the distribution of singularity exponents. We provided a detailed evaluation of our algorithm and compared it against several alternative approaches using the full Train and Test sets of TIMIT. The comparisons showed that our algorithm is more accurate than state-of-the-art ones. These encouraging results with a very simple analysis of singularity exponents confirm that nonlinear dynamics of the speech signal can indeed be geometrically accessed through the use of the MMF.

Fig. 1 -

 1 middle displays the time-evolution of local histograms of h[n]. Each column in this figure shows a histogram of SE values computed on a window of 32ms length around the time t. The red color shows the highest probability and the blue color shows zero probability. One can easily observe that there is a change in the position of peak and in the variability of the distribution at phoneme boundaries.

Figure 1 :

 1 Figure 1: (top) A part of a speech signal from TIMIT database. Manually-positioned phoneme boundaries are marked with vertical red lines. (middle) The local histogram of singularity exponents in a 32ms window around each time instant. Red corresponds to maximum probability and dark blue corresponds to zero probability. (bottom) The number of points belonging to the MSM in a sliding window of 32ms length.

Figure 2 :

 2 Figure 2: The speech signal "Masquerade parties tax one's imagination" from the TIMIT database and the ACC functional of Eq. (7).

Figure 3 :

 3 Figure 3: (top) Two examples of ACC for the original signal and (bottom) The ACC for the low-passed filtered signal. Phoneme boundaries are marked with vertical red lines.

Figure 4 :

 4 Figure 4: The sensitivity of the SE-LLRT to threshold.

  the MSE of approximating ACC by a straight line from ACC[m i ] to ACC[m j ].

	12: 13: end while k 1 ← k 1 + 1;
	14:
	15:	⊲ E	m i →m j 1

is 16:

Table 2 :

 2 The comparative table of segmentation results. The scores are reported as percentages.

	tolerance score Dusan et al SE-ACC-Train SE-LLRT-Train SE-LLRT-Test
		HR	22.8	31.7	32.3	32.4
	5ms	FA	79.7	70.2	68.5	62.02
		OS	12.8	6.4	2.5	4.61
		HR	-	52.8	53.6	53.16
	10ms	FA	-	50.4	47.7	49.18
		OS	-	6.4	2.5	4.61
		HR	59.2	65.5	66.3	65.39
	15ms	FA	47.5	38.4	35.4	37.50
		OS	12.8	6.4	2.5	4.61
		HR	-	72.4	72.5	71.7
	20ms	FA	-	31.94	29.2	31.5
		OS	-	6.42	2.5	4.61
		HR	75.3	76.2	76.1	75.17
	25ms	FA	33.2	28.3	25.8	28.14
		OS	12.8	6.4	2.5	4.61
		HR	-	78.8	80.5	77.41
	30ms	FA	-	26	23.7	26.00
		OS	-	6.4	2.5	4.61
	Table

Table 3 :

 3 The comparative table of global performance measures.

	tolerance	score	Dusan et al SE-ACC-Train SE-LLRT-Train SE-LLRT-Test
	5ms	R-value F 1 -value	0.29 0.21	0.39 0.31	0.41 0.32	0.41 0.32
	10ms	R-value F 1 -value	--	0.57 0.51	0.60 0.53	0.58 0.52
	15ms	R-value F 1 -value	0.60 0.55	0.68 0.63	0.70 0.65	0.69 0.64
	20ms	R-value F 1 -value	--	0.74 0.70	0.76 0.72	0.74 0.70
	25ms	R-value F 1 -value	0.73 0.71	0.77 0.74	0.79 0.75	0.77 0.735
	30ms	R-value F 1 -value	--	0.79 0.76	0.81 0.77	0.79 0.76

Table 4 :

 4 The segmentation results of[START_REF] Esposito | Text independent methods for speech segmentation[END_REF], for three different coding schemes. Tolerance : 20ms.

	Coding scheme	(a,b,c)	HR	OS	FA	R-value F-value
	8-Mel-bank	(2,optimal value,5) 82 18.32 30.69	0.743	0.75
	5-MFCC		76 12.31 31.33	0.736	0.72
	Log Area Ratio		70	6.31 34.16	0.72	0.68

with the two lines that connects ACC[START_REF] Kaiser | Some observations on vocal tract operation from a fluid flow point of view[END_REF], ACC[c 1 ] and ACC[n 1 ]

Acknowledgement

The first author is funded by the INRIA CORDIS doctoral program.


Arthur Hellouin De Menibus
email: arthur.hellouin-de-menibus@cea.fr

Quentin Auzoux

Jacques Besson

Jérôme Crépin

Temperature Increase of Zircaloy-4 Cladding Tubes due to Plastic Heat Dissipation during Tensile Tests at 0.1 to 10 s -1 Strain Rates

Keywords: Zirconium, Zircaloy-4, Plastic Dissipation, Infrared, Reactivity Initiated Accident

This study is focused on the impact of rapid Reactivity Initiated Accident (RIA) representative strain rates (about 1 s -1 [1]) on the behavior and fracture of unirradiated Cold Work Stress Relieved Zircaloy-4 cladding tubes. Uniaxial ring tests (HT) and plane strain ring tensile tests (PST) were performed in the 0.1 to 10 s -1 strain rate range, at 25 °C. The local temperature increase due to plastic dissipation was measured with a high-speed infrared camera.

Limited temperature increases were measured at 0.1 s -1 strain rate. Limited but not strongly localized temperature increases were measured at 1 s -1 . Large temperature increase were measured at 5 and 10 s -1 (142 °C at 5 s -1 strain rate in HT tests). The local temperature increase induced heterogeneous temperature fields, which enhanced strain localization and resulted in a reduction of the plastic elongation at fracture.

INTRODUCTION

The hypothetical control rod ejection in a Pressurized Water Reactor leads to a rapid RIA, which triggers rapid thermal expansion of the fuel pellets. These pellets expansion induced a strain-controlled loading on the cladding at approximately 1 s -1 circumferential strain rate [START_REF]Nuclear fuel behaviour under reactivity-initiated accident (RIA) conditions[END_REF]. This phase is referred to as PCMI (Pellet Cladding Mechanical Interaction). Such a 1 s -1 strain rate is susceptible to cause an important increase of the cladding temperature due to plastic heat dissipation. The aim of the present study is to quantify this temperature increase and its effect on the Zircaloy-4 cladding tubes behavior and fracture.

The temperature increase due to plastic heat dissipation assuming adiabatic conditions can be estimated for uniaxial ring tensile tests (Figure 1) performed on unirradiated Zircaloy-4, using the following equation: ̇ ̇ (eq. 1)

where is the uniaxial stress (MPa), ̇ is the plastic strain rate (s -1 ), is the material density (kg/m 3 ), is the specific heat (J/kg/K) and is the Taylor-Quinney coefficient [START_REF] Taylor | The latent energy remaining in a metal after cold working[END_REF][START_REF] Quinney | The emission of the latent energy due to previous cold working when a metal is heated[END_REF]. The β coefficient varies with plastic strain [START_REF] Macdougall | Determination of the Plastic Work Converted to Heat Using Radiometry[END_REF] and can be greater than one. Indeed, the heat released during a given time period is the plastic energy dissipated into heat plus a possible release of energy previously stored in the material [START_REF] Vivier | On the stored and dissipated energies in heterogeneous rate-independent systems[END_REF]. Regarding pure zirconium, it was shown that the β coefficient depended on the material crystallographic texture: it is higher in non-basal than in basal directions [START_REF] Padilla | Spatiotemporal Thermal Inhomogeneities During Compression of Highly Textured Zirconium[END_REF]. More precisely, the results in [START_REF] Padilla | Spatiotemporal Thermal Inhomogeneities During Compression of Highly Textured Zirconium[END_REF] suggest that setting β=90% is an acceptable assumption when stretching the samples in directions parallels to zirconium <a> axis. In directions parallels to <c> axis, β increases from 50% to 75% from 0% to 10% plastic strain, and remains almost constant for higher plastic strain. In CWSR Zircaloy-4 cladding tubes, the crystallographic texture is such that <c> axis are mostly in the radialcircumferential plane at 30° to the radial direction [START_REF] Murty | Texture development and anisotropic deformation of Zircaloys[END_REF]. Therefore, the axial direction corresponds to a non-basal direction, following which β=90%. The circumferential direction is not purely non basal, one may assume it results in β=50%-90%. To estimate an upper bound for the temperature increase, a value of β=90% was used and the stress level was assumed to be constant and equal to the maximum stress. Therefore, equation 1 can be integrated from zero to a given plastic strain level ( p ) considering β is independent of plastic strain:

(eq. 2)

The temperature increase at fracture is about 270 °C using the maximum stress and the plastic elongation at fracture:

(eq. 3)
With the elongation at fracture, the calibrated gauge length equal to 3 mm, the stress at fracture and E the elasticity modulus.

Considering the maximum duration of a CABRI test (76.4ms for REP-Na4 [START_REF] Papin | {Summary and interpretation of the CABRI REP-Na program}[END_REF]), the thermal diffusion distance √ , with the thermal diffusivity of zirconium (m 2 /s) [START_REF]Thermophysical properties of materials for nuclear engineering: a tutorial and collection of data[END_REF] and the time (s), is about 0.7 mm. Thus, the diffusion distance is similar to the sample dimensions (0.57 mm in thickness for instance). As a consequence, adiabatic conditions are not fulfilled in RIA conditions. The local maximum temperature increase due to plastic heat dissipation is expected to decrease when the pulse duration increases. The effect of a local temperature increase on the fracture strain (at local and at macroscopic scale) is not trivial. On one hand the temperature increase may increase the material ductility, but on the other hand the induced heterogeneous thermal field could promote strain localization.

In order to measure the local heating at RIA strain rates and its potential impact on the behavior and fracture, uniaxial tensile test on ring and plane strain tensile tests were performed at strain rates ranging from 0.1 to 10 s -1 at 25 °C on Zircaloy-4 cladding tubes and were monitored with a high-speed infrared camera. First, the thermal fields measured during such tests are described and their impact on the mechanical behavior and fracture is discussed. Then, the temperature increase at higher temperature (350 °C and 480 °C) is theoretically estimated using the energy balance.

MATERIAL

The studied material consisted in 9.5 mm external diameter and 0.57 mm thickness tubes of Cold Worked Stress Relieved (CWSR) unirradiated and non hydrided Zircaloy-4 provided by CEZUS Company ® . The grain size measured by statistical analysis of polarized light optical pictures was 2.5±0.7 µm in radial, 3.7±0.7 µm in circumferential and of about 10 µm in the tube axial direction. This material was obtained with the standard manufacturing route for CWSR Zircaloy-4. For this route, most of the zirconium <c> axis are oriented at ±30° to the tube radial direction and one of the ̅ direction is aligned with the tube axial direction [START_REF] Murty | Texture development and anisotropic deformation of Zircaloys[END_REF].

PROCEDURES a. Tensile tests

Uniaxial ring (Figure 1) and plane strain ring tensile tests (Figure 2) were performed. In the following, they would be referred to as HT for Hoop Tensile test and PST for Plane Strain Tensile tests, respectively. The HT gauge section was 2 mm wide (W) and 3 mm (L 0 ) long. The PST geometry has no clearly calibrated gauge length because stress and strain fields are heterogeneous in PST samples. However, the gauge section was considered equal to the notches diameter L 0 =2 mm in order to derive from the displacement a macroscopic strain value. During a PST test at 25 °C, it was shown by Le Saux [START_REF] Saux | Comportement et rupture de gaines en Zircaloy-4 détendu vierges, hydrurées ou irradiées en situation accidentelle de type RIA[END_REF] that a centered area of 1 mm in circumferential direction and 2 mm in axial direction reaches plane strain conditions. Both types of tensile tests were done using two 8.30 mm diameter D-shape mandrels on an INSTRON servo-hydraulic tensile machine (Figure 3). No "dog-bone" part was used and no lubricant was inserted between mandrel and ring. Above 1 kN, the accuracy of the load measurements is better than 0.1%, and the discrepancy on the imposed strain rate is lower than 0.3%. The load-displacement curve was recorded with a 1024 Hz frequency.

b. Infrared monitoring

The tests were monitored with a FLIR SC7000 high-speed infrared camera at frame rates up to 3600 Hz. The spatial resolution of the infrared measurement was 70 µm/pixel in HT tests and 110 µm/pixel in PST tests. Before the tensile tests, a 7778 BBQ black RUSTOLEUM ® black paint coating was hand sprayed in a reproductive manner. The coating thickness was ranging from 10 µm to 15 µm based on scanning electron microscope cross sections observations. Comparison of infrared camera and thermocouple measurements showed that the coating emissivity was at room temperature and was constant up to 200 °C. It was checked that a 5% uncertainty on the emissivity results in a 0.4 °C uncertainty on the absolute temperature at 25 °C, and a 1.9 °C uncertainty at 90 °C. The HT tests were monitored with the infrared camera at a 30° angle from the normal to the sample gauge section, and PST tests with a 0° angle. The reason of the tilted angle is HT samples were also monitored with a high speed optical camera for DIC application. The tilted angle is not an issue regarding the temperature measurements because it was shown by Poncelet [11] that the emissivity of highly emissive black coatings remained constant up to a 45° angle. In the following, all the thermal measurements are expressed as relative temperature ( ) to get rid of any parasitic reflections that were supposed not to vary during the tests short durations (the longest test lasted less than 5.5 s).

c. Thickness reduction at fracture and elongation at fracture measurements

The thickness reduction at fracture Δe/e was measured with an optical macroscope on the radial axial plane, with 20 measurements on both sides of the fractured specimens. The uncertainty on the average plastic strain in the thickness direction was estimated to be 1.5% for PST samples and 0.6% for HT samples. The HT thickness reduction was measured in the center of the 2 mm gauge section, in a 0.3 mm axially extended area.

Plastic elongation at fracture was calculated from the load-displacement curves and the gauge lengths defined above.

RESULTS

HT tests and PST tests both require a limited quantity of material. The loading biaxiality level in HT tests is mainly uniaxial in the circumferential direction, which not representative of the loading biaxiality expected in RIA conditions that is from plane strain to equibiaxial conditions. Thus, PST tests, which allow reaching plane strain conditions, were also performed. The mechanical test results and associated temperature measurements are listed for each individual test in Table 1. Averaged values for each type of tests are presented in Table 2. The parameter is the yield stress (MPa), is the maximum stress (MPa), (%) is the plastic strain at the maximum load, is the thickness reduction at fracture, is the average temperature increase at fracture on a gauge section centered area (1.4 x 3.2 mm 2 for HT tests and 2.0 x 0.9 mm 2 for PST tests) and is the maximum temperature increase at fracture, measured on the hottest pixel on the whole infrared frame. a. HT tests at 25 °C During HT tests at 25 °C, the temperature increase visible on the infrared frames was less localized at 0.1 s -1 than at 5 s -1 due to heat conduction (Figure 4). The infrared detector saturated at 125 °C (100 °C temperature increase). When the temperature increases at fracture were greater than 100 °C, for instance in HT test at 5 s -1 strain rate, the temperature increase evolution was extrapolated following the linear trend of the temperature increase observed since the maximum load. The maximum temperature increase at fracture reached 142 °C at 5 s -1 strain rate, while it was limited to 12 °C at 0.1 s -1 . This high local temperature enhanced the material ductility: the thickness reduction at fracture increases slightly when the strain rate increased from 0.1 to 5 s -1 (Table 2). On the other hand, the heterogeneous temperature field enhanced strain localization in the gauge section: the plastic elongation at fracture ( ) decreased from 25% to 18% (-29%). The strain in the gauge section is clearly more localized at 5 s -1 than at 0.1 s -1 macroscopic strain rate (Figure 5). In addition, higher strain rates led to a steeper load decrease right after the load maximum (Figure 6). A zoom on the elastic phase of the HT tests showed a thermoelastic effect (Figure 7). The specimen temperature decreased during the tensile loading in the elastic regime. Under adiabatic conditions, the theoretical temperature decrease can be calculated with the following equation:

(eq. 4)

with α the expansion coefficient (K -1 ), calculated as a weighting of the expansion in zirconium < a > and <c> directions [START_REF]Thermophysical properties database of materials for light water reactors and heavy water reactors[END_REF]:

(eq. 5) (eq. 6) with the typical Kearns factors [START_REF] Murty | Texture development and anisotropic deformation of Zircaloys[END_REF] for CWSR claddings equal to in the radial direction and in the circumferential one. The term T is the average temperature (room temperature in the present case, in K), E is the elasticity modulus (GPa), ρ is the material density (kg/m 3 ), C is the specific heat (J/kg/K), Δε is the elastic strain variation. As expected, the measured thermoelastic effect was lower at 0.1 s -1 than at 5 s -1 , which was itself lower than the adiabatic estimate due to heat losses by conduction (the elasticity contribution of the tensile machine and setup was assumed to be negligible). The discrepancy between experimental measurements and the theoretical estimate might also be due to the non-uniform stress distribution in the specimen due to some bending up to 0.5% plastic strain [START_REF] Saux | Comportement et rupture de gaines en Zircaloy-4 détendu vierges, hydrurées ou irradiées en situation accidentelle de type RIA[END_REF]. The decrease in temperature at 5 s -1 was limited to -0.4 °C, which was negligible in comparison to the temperature increase due to plasticity. The good reproducibility of the measurement of the thermoelastic induced temperature variation demonstrated that the setup, procedures and temporal synchronization are sufficiently accurate to capture such a limited effect in a reproducible fashion. b. PST tests at 25 °C During PST tests at 25 °C, similarly to HT tests, the temperature increase was less pronounced at 0.1 s -1 strain rates than at 5 s -1 (Figure 8). The temperature distribution measured at 5 s -1 was in agreement with the circumferential strain map measured by digital image correlation by Le Saux [START_REF] Saux | Comportement et rupture de gaines en Zircaloy-4 détendu vierges, hydrurées ou irradiées en situation accidentelle de type RIA[END_REF]: the circumferential strain is localized at the notches and in two bands that surrounded the central area (the central area is in plane strain).

The maximum temperature increase at fracture was 75 °C at 5 s -1 macroscopic strain rate, while it was about 5 °C at 0.1 s -1 strain rate (Figure 9). Some specimens failed on the opposite side to the one monitored by the infrared camera, the PST sample tested at 5 s -1 showed in Figure 8(b) for instance. The temperature increase was similar whatever the side on which the sample was fractured, as shown by the markers "front" and "back" in Figure 9, which identify the samples that failed in front of the camera or on the opposite side. The temperature increase was mainly located at the notches, where the stress biaxiality was low (uniaxial) while the central area in plane strain conditions showed limited temperature increase.

The fracture of the PST samples tested in the present study were consistent with the scenario described by Le Saux [START_REF] Saux | Comportement et rupture de gaines en Zircaloy-4 détendu vierges, hydrurées ou irradiées en situation accidentelle de type RIA[END_REF]. First, a through thickness plastic strain instability initiated on the inner side in the center of the gauge section and propagated in a plane at 45° to the radial direction. Secondly, the crack propagated axially toward the notches and bifurcated before the complete fracture of the sample. It was observed in the present study that the aspect of the fractured samples was identical whatever the strain rate. Thus, the fracture of PST sample was not initiated where the temperature increase due to plastic dissipation was the highest (at the notches), but were the stress triaxiality level was the highest (in the plane strain area). Consequently, the thickness reduction at fracture was independent of the strain rate. Similar observation was reported in [START_REF] Link | Failure of Zircaloy cladding under transverse plane-strain deformation[END_REF] on PST samples when the strain rate was increased from 10 -3 s -1 to 10 2 s -1 . On the contrary, the plastic elongation at fracture was affected by the temperature increase at the notches. It decreased in -20% when the strain rate was increased from 0.1 to 5 s -1 .

As observed in HT tests, a higher strain rate led to a steeper load decrease after the maximum load (Figure 10). The magnitude of the slope was only 1.5 times steeper when the strain rate was increased from 0.1 to 5 s -1 while it was 2.1 times steeper in HT tests. The load maximum was 974 MPa at 0.1 s -1 , increased up to 1020 MPa at 1 s -1 , and increased up to 1029 MPa at 10 s -1 . The temperature distributions measured between the two notches for PST specimens tested at different strain rates, after the load maximum are reported in Figure 11. For the PST specimen tested at 0.1 s -1 , the temperature distributions were homogeneous and showed almost no local heating. At 1 s -1 , the temperature increased and heat conduction resulted in a relatively homogeneous distribution. When the strain rate was further increased to 2.5, 5 and 10 s -1 , the thermal diffusion length decreased and the temperature increased sharply at the notches where high plastic strains were reached. Therefore, the temperature on the gauge section center first increased and then declined when the strain rate increased from 0.1 to 10 s -1 . This is also visible in the evolution of the average temperature increase in Figure 9(a). The studied strain rate range turned out to be a double transition from almost no local heating to important temperature increase and from significant heat conduction to nearly adiabatic conditions. c. HT at higher temperature: 350 °C and 480 °C Infrared measurements were also performed at higher temperatures, but the black coating cracked before reaching the maximum load. Thus, the heating due to plasticity dissipation was estimated theoretically at 350 °C and 480 °C at 5 s -1 strain rate in HT tests. The procedure consisted in computing the ratio of the density K in kg/m 3 [START_REF]Thermophysical properties of materials for nuclear engineering: a tutorial and collection of data[END_REF] times the specific heat [START_REF]Thermophysical properties database of materials for light water reactors and heavy water reactors[END_REF] in J/kg/K and times the average temperature increase divided by the integral of the macroscopic stress times the plastic strain ∫ , as shown in the following equation:

∫ ∫ ∫ (eq. 7)
The subscript 25, 350 and 480 correspond to 25 °C, 350 °C and 480°C, respectively. Additional tensile tests were performed at 350 °C and 480 °C, and then were used to determine the evolution of ∫ and ∫ in function of the plastic elongation. The only unknown parameters were the temperature increases at high temperature ( and ). The equation 7 was deduced from the two following assumptions: the β coefficient remained constant with the strain level and the temperature, and heat losses are small compared to the heat production. The estimate obtained using equation 7, showed in Figure 12, gave an average temperature increase at fracture equal to and , while it was during room temperature testing at a 17.5% plastic elongation (close to the strain at fracture at 25 °C). This lower average temperature increase during higher temperature tests was mainly due to the lower stress level for a given plastic strain.

DISCUSSION: CONSEQUENCES OF THE LOCAL TEMPERATURE INCREASE ON THE ANALYSIS OF MECHANICAL TESTS AT RIA-PCMI PHASE REPRESENTATIVE STRAIN RATES

Le Saux et al. [START_REF] Saux | A model to describe the anisotropic viscoplastic mechanical behavior of fresh and irradiated Zircaloy-4 fuel claddings under RIA loading conditions[END_REF] analysis of the PROMETRA database [START_REF] Cazalis | The PROMETRA program: Fuel cladding mechanical behavior under high strain rate[END_REF] showed that the stress strain sensitivity ̇ , represents the material capability to resist to strain localization, increased from 0.013 at 25 °C to 0.022 at 400 °C and then rose steadily to about 0.24 at 800 °C. In addition, the thermal conductivity increases in 30% from 25 °C to 480 °C for Zircaloy alloys [START_REF]Thermophysical properties database of materials for light water reactors and heavy water reactors[END_REF], thus limiting the localization of the temperature increase. These two last statements together with the lower temperature increase at high temperature indicate that the magnitude of the elongation at fracture reduction due to material autoheating would likely be lower at 350 °C than at 25 °C, and even lower at 480 °C.

The cladding circumferential strain rate in RIA conditions is about 1 s -1 . Therefore, limited local heating (<25 °C) due to plastic dissipation may occur after the maximum load and the increase in temperature should not be strongly localized due to significant heat conduction. Nonetheless, a limited localized temperature increase could enhance strain localization. Now that the general trend and magnitude of the temperature increase is known in the strain rate range from 0.1 to 10 s -1 , further studies should focus on strain rates close to 1 s -1 , at RIA representative temperatures. Due to heat conduction, an accurate estimate of the local temperature cannot be done assuming adiabatic conditions and would require a coupled thermal-mechanical calculation scheme.

Concerning the interpretation of previous experimental results, the present study showed the local temperature increase becomes significant after the maximum load is reached. Thus, the parameters such as the maximum load and uniform strain would not be significantly affected. On the other hand, the plastic elongation at fracture should not be used to quantify the material ductility at strain rate greater than 0.1 s -1 because of the local temperature increase induced by the plastic heat dissipation. Indeed, the -29% reduction in the plastic elongation at fracture measured during HT tests at 0.1 and 5 s -1 was not due to a material ductility sensitivity to strain rate. It was due to the strain localization induced by the heterogeneous thermal field, which depends on the geometry of the tensile test sample.

This last point should be reminded of, when the elongation at fracture obtained by uniaxial ring testing is converted to higher stress biaxialities using correction factors [START_REF] Jernkvist | A strain based clad failure criterion for reactivity initiated accidents in light water reactors[END_REF]. The present study showed that this method is not directly applicable for strain rate greater than 0.1 s -1 (and a fortiori for strain rate higher than 1 s -1 ). A more local measurement of the fracture strain, such as the thickness reduction or full field techniques such as Digital Image Correlation should be used instead of the plastic elongation at fracture.

CONCLUSION

Infrared measurements of the local heating due to plastic dissipation were performed during tensile tests at 25 °C in the strain rate range from 0.1 to 10 s -1 with two different tensile tests on CWSR Zircaloy-4 cladding tubes. Limited temperature increase was measured at 0.1 s -1 strain rate, while it reached 142 °C at 5 s -1 strain rate in uniaxial ring tests. Local temperature increase was observed between 0.1 and 1 s -1 . Nevertheless, the temperature increase was limited in magnitude (< 25°C in PST tests at 1 s -1 at fracture) and was not strongly localized due to thermal conduction. At strain greater than 1 s -1 , the local temperature increase induced heterogeneous temperature fields. It induced a plastic elongation at fracture reduction. The local thickness reduction measurement showed this was not due to a drop in material ductility. Thus, fracture strains of mechanical tests at strain rates equal or greater than 1 s -1 should be measured at a local scale to assess the material intrinsic fracture behavior.

Figure 1 :

 1 Figure 1: HT sample drawing.

Figure 2 :

 2 Figure 2: PST sample drawing.

Figure 3 :

 3 Figure 3: Black painted PST sample placed on the D-shape mandrels.

Figure 4 :

 4 Load-displacement evolution and infrared frames in a HT test at (a) 0.1 s -1 and (b) 5 s -1 strain rate. Several temperature scales were used.

Figure 5 :

 5 Figure 5 : Aspect of additional HT samples gauge length tested at 25 °C on the last frame prior to fracture (6250Hz). Prior to testing, samples surfaces were sputtered with black and white paint drops for digital image correlation purpose (not detailed in the present paper).

Figure 6 :

 6 Figure 6: Individual load-displacement in HT tests with the (a) average temperature and (b) maximum temperature evolution, and resulting average evolution.

Figure 7 :

 7 Figure 7: Zoom on the elastic phase showing a thermoelastic phenomenon -load-displacement, maximum and average temperature evolutions during a HT test at 25 °C and 0.1 s -1 and 5 s -1 strain rates.

Figure 8 :

 8 Load-displacement evolution and infrared frames during a PST test at (a) 0.1 s -1 and (b) 5 s -1 strain rate. Several temperature scales were used.

Figure 9 :

 9 Figure 9: Individual load-displacement in PST tests with the (a) average temperature and (b) maximum temperature evolution, and resulting average evolution.

Figure 10 :

 10 Figure 10 : Zoom on the maximum load of individual load-displacement evolution in PST tests and resulting average evolution. The first arrow shows the load increase with strain rate. The second arrow shows the increase in slope of the load-plastic strain evolution after the load maximum when the strain rate is increased.

Figure 11 :

 11 Figure 11: Temperature at 0.32 plastic elongation, soon after the maximum load, along a line joining the two notches of a plane strain specimen.

Figure 12 :

 12 Figure 12: Experimental Load-plastic strain curves of HT tests at 25°C, 350°C and 480°C, experimental average temperature increase (ΔT avg )-plastic strain curve at 25°C and calculated at 350 and 480°C.

Table 1 : Mechanical results and infrared camera measurements of the individual tensile tests performed (initial test temperature is 25 °C).

 1

					Mechanical results			Infrared camera measurements
		Test type	Strain rate (s -1 )	σy0.2 (MPa)	σu (MPa)	(%)	(%)	(%)	At maximum load	(°C) fracture At	At maximum load	(°C) At Fracture
		HT	0.1	778	855	2.7	24.7	34.5	1.5	6.2	1.3	11.9
		HT	0.1	757	846	3.1	15.2	33.3	1.7	6.3	1.8	12.8
		HT	5	839	903	1.5	16.5	36.9	1.3	52.9	2.7	142.0
		HT	5	842	911	1.6	18.2	35.4	0.1	48.1	0.7	130.8
		HT	5	825	902	1.9	19.1	36.0	0.3	57.5	2	140.9
		PST	0.1	808	981	4.2	15.2	13.4	1.5	3.8	1.5	5.9
		PST	0.1	788	981	5.1	16.8	20.1	1.7	4.3	1.8	4.8
		PST	1	829	993	4.4	14.0	13.8	4.7	22.5	12.0	30.3
		PST	1	842	1009	4.6	15.3	14.8	4.9	23.4	10.4	29.1
		PST	2.5	848	1006	4.5	11.3	15.5	5.4	29.4	21.0	46.1
		PST	5	854	1030	4.3	12.5	17.7	1.9	29.6	17.2	59.8
		PST	5	842	1031	3.9	14.1	17.2	1.9	46.3	15.6	74.2
		PST	5	848	1024	4.5	12.5	13.3	1.6	25.2	19.3	62.0
		PST	5	791	1025	3.9	14.0	20.6	1.9	39.5	14.4	67.1
		PST	5	832	1014	4.2	13.1	15.2	1.3	41.0	15.5	67.9
		PST	5	850	1018	4.2	10.9	17.1	1.3	22.7	13.4	28.6
		PST	10	869	1036	3.5	10.9	15.0	2.4	28.2	25.7	73.6
		PST	10	806	1032	3.5	12.5	20.4	0.9	37.5	21.8	84.7
					Mechanical results				Infrared camera measurements
	Test type	Strain rate (s -1 )	Number of tests	(MPa) σy0.2	σu (MPa)	(%)	(%)	(%)		(°C)		(°C)

At maximum load At fracture At maximum load At Fracture HT 0.1 2 768±15 851±6 2.9±0.3 25.0±0.4 33.9±0.8 1.6±0.1 6.3±0.1 1.6±0.4 12.4±0.6 HT 5 3 835±9 905±5 1.7±0.2 17.9±1.3 36.1±0.8 0.6±0.6 52.8±4.7 1.8±1.0 137.9±6.2 PST 0.1 2 798±14 981±0 4.7±0.6 16.0±1.1 16.8±4.7 1.6±0.1 4.1±0.4 1.7±0.2 5.4±0.8

	PST	1	2	836±9	1001±11	4.5±0.1	14.7±0.9 14.3±0.7	4.8±0.1	23.0±0.6	11.2±1.1	29.7±0.8
	PST	2	1	848	1006	4.5	11.3	15.5	5.4	29.4	21.0	46.1
	PST	5	6	836±23	1024±7	4.2±0.2	12.9±1.2 16.9±2.5	1.7±0.3	34.1±9.5	15.9±2.1	63.3±8.8
	PST	10	2	838±45	1034±3	3.5±0.0	11.7±1.1 17.7±3.8	1.7±1.1	32.9±6.6	23.8±2.8	79.2±7.8

Table 2 : Mechanical tests and infrared camera measurements averaged results.

 2

ACKNOWLEDGMENTS

The authors are grateful to Abdennour Meddour and Jean-Christophe Teissedre from MinesParistech for technical assistance and Herve Louche from Université Montpellier 2 for enlightening advices.


K Motchon
email: kmd.motchon@ed.univ-lille1.fr

K Pekpe

J-P Cassar

S De Bièvre

On the Distinguishability of Positive Linear Time-Invariant Systems with Affine Parametric Uncertainties

Keywords: Distinguishability, Uncertain system, Parametric uncertainty, Positive system

Studies of distinguishability have focused on the case of Linear Time-Invariant systems without uncertainties. In this work, distinguishability is studied for Positive Linear Time-Invariant systems with affine parametric uncertainties in the state space model. We propose a definition of distinguishability adapted to this new context and give a characterization of this notion. The approach used is based on the estimate of the reachable output space of the systems. Under suitable assumptions, a sufficient condition for distinguishability is established.

INTRODUCTION

Consider the two uncertain Linear-Time Invariant (LTI) systems described by the following set of differentialalgebraic equations :

S k ẋk (t) = A k (θ k ) x k (t) + B k u(t), y k (t) = C k x k (t), x k (0) = x o k ∈ X o k , (1)
with k = 1, 2 and where x k (t) ∈ E k ⊆ R n , y k (t) ∈ R m and u(t) ∈ R l are respectively the state vector, the output vector and the input vector of the system S k ; E k is the state space of S k and the set X o k is such that X o k ⊆ E k ; the matrix B k is the input matrix and C k is the output matrix ;

θ k = [θ k,1 θ k,2 • • • θ k,p k ]
T is a vector of uncertain real parameters θ k,r , r = 1, 2, . . . , p k and the state matrix A k (θ k ) has the form

A k (θ k ) = A k + θ k,1 A k,1 + θ k,2 A k,2 + • • • + θ k,p k A k,p k ,
where A k , A k,1 , . . . , A k,p k are known matrices. The lower bound θ k,r and the upper bound θ k,r of each real parameter θ k,r ∈ [θ k,r ; θ k,r ] are assumed to be known. The systems S 1 and S 2 can represent two distinct operating modes of an affine switched system with uncertain parameters. This paper is concerned with the property of distinguishability between the uncertain systems S 1 and S 2 . Most distinguishability problems are discussed for LTI systems without uncertainties [START_REF] Avdeenko | The problem of distinguishability of state space models[END_REF], [START_REF] Cocquempot | Fault detection and isolation for hybrid systems using structured parity residuals[END_REF], [START_REF] Grewal | Identifiability of linear and nonlinear dynamical system[END_REF], [START_REF] Lou | The distinguishability of linear control systems[END_REF], [START_REF] Lou | Conditions for distinguishability and observability of switched linear systems[END_REF], [START_REF] Motchon | Operating modes distinguishability condition in switching systems[END_REF] S. De Bièvre is supported in part by the Labex CEMPI (ANR-11-LABX-0007-01). The work of K. Motchon is supported by a PhD Student Scholarship from Lille 1 University.

and such studies have applications in Fault Detection and Isolation (FDI) [START_REF] Bayoudh | Coupling continuous and discrete event system techniques for hybrid system diagnosability analysis[END_REF], [START_REF] Cocquempot | Fault detection and isolation for hybrid systems using structured parity residuals[END_REF]), and in active mode detection in switched affine systems [START_REF] Domlan | Active mode estimation for switching systems[END_REF], [START_REF] Hakem | Parameter-free method for switching time estimation and current mode recognition[END_REF]). But it is well known that to describe a system with a linear model, approximations have generally to be made leading to uncertainties in the parameters of the linear model. The motivation of this study is to take such uncertainties into account in the analysis of the distinguishability property. In the sequel,

• U is a sub-vector space of L 1 ([0 ; T ], R m ). • U ⊆ L 1 ([0 ; T ], R m )
denotes the space of admissible inputs of the systems S 1 and S 2 .

• X o = X o 1 × X o 2 denotes the set of all admissible initial state pairs (x o 1 , x o 2 ) of systems S 1 and S 2 . • [θ k ; θ k ] := [θ k,1 ; θ k,1 ] × • • • × [θ k,p k ; θ k,p k ] denotes
the set of all admissible values of the vector θ k .

• Θ := [θ 1 ; θ 1 ] × [θ 2 ; θ 2 ] ⊂ R p1 × R p2 is the set of all
admissible parameter pairs (θ 1 , θ 2 ) of systems S 1 and S 2 .

The output vector y k (t) of S k can be expressed as follows (Larminat ( 2007)):

y k (t) = C k e t A k (θ k ) x o k + t 0 e (t-τ ) A k (θ k ) B k u(τ ) dτ .
Consequently, y k (t) is a function of the initial state vector x o k , the input signal u and the vector θ k . This function is denoted by y k (t) ≡ y k (t, x o k , u, θ k ). For a given value θ k ∈ [θ k ; θ k ] of the vector θ k , let S(θ k ) be the system defined as follows:

S(θ k ) : ẋk (t) = A k (θ k ) x k (t) + B k u(t), y k (t) = C k x k (t), x k (0) = x o k .
According to the definition proposed in [START_REF] Lou | The distinguishability of linear control systems[END_REF], two systems S(θ 1 ) and S(θ 2 ) are distinguishable on [0 ; T ] if and only if for all triplets (

x o 1 , x o 2 , u) ∈ R n × R n × U such that (x o 1 , x o 2 , u) = (0, 0, 0), the outputs y 1 (•, x o 1 , u, θ 1 ) and y 2 (•, x o 2 , u, θ 2 ) are not identical on [0 ; T ].
In particular, given u = 0, no output of S(θ 1 ) can be equal to an output of S(θ 2 ). A generalization of this notion to the case of LTI systems with affine uncertainties is the (X o , U, Θ)distinguishability that we define as follows: Definition 1. Systems S 1 and S 2 are said to be (

X o , U, Θ)- distinguishable on [0 ; T ] if for all (x o 1 , x o 2 ) ∈ X o , all u ∈ U and all (θ 1 , θ 2 ) ∈ Θ with (x o 1 , x o 2 , u) = (0, 0, 0), the functions y 1 (•, x o 1 , u, θ 1 ) and y 2 (•, x o 2 , u, θ 2 ) are not identical on [0 ; T ]. If not, the two systems are said to be (X o , U, Θ)- indistinguishable on [0 ; T ].
This paper proposes a sufficient condition for (X o , U, Θ)distinguishability. The condition is based on an estimate of the reachable output spaces of S 1 and S 2 which is obtained with a method developed in [START_REF] Kieffer | Guaranteed nonlinear state estimation for continuous-time dynamical models from discrete-time measurements[END_REF] and [START_REF] Meslem | Atteignabilité hybride des systèmes dynamiques continus par analyse par intervalles : Application à l'estimation ensembliste[END_REF]. When S 1 and S 2 are positive systems and under some assumptions on the spaces X o k , U and the matrices A k,r , it is straightforward to obtain this estimate.

In Section 2 some basic notations and definitions are given and the estimation technique is recalled. Section 3 is devoted to the proof of the sufficient condition for distinguishability. An illustration of this result with an academic example is given. Some remarks and open problems conclude the work.

PRELIMINARIES

Some basic notation and definitions

Let A = (a ij ) ∈ R n×m be a matrix and, v, w, w and w be four vectors of R p .

• A is said to be a Metzler-matrix if ∀i = j, a ij ≥ 0.

• A is said to be a non-negative matrix (A 0) if ∀i, j, a ij ≥ 0 and R n×m + := {M ∈ R n×m : M 0} denotes the set of all non-negative matrices of R n×m .

• v is said to be a non-negative vector (v 0) if ∀i,

v i ≥ 0 and R p + := {ξ ∈ R p : ξ 0} denotes the set of all non-negative vectors of R p . • v w if w -v 0 ; w ∈ [w ; w] if w w and w w. • v 0 if v 0 and ∃i ∈ {1, 2, . . . p} such that v i > 0. • v w if v -w 0.

Positive linear dynamical systems : definition and characterization

Positive dynamical systems are a class of dynamical systems which is studied in detail in [START_REF] Farina | Positive Linear Systems: Theory and Applications[END_REF]. We only recall the definition and one basic property of these systems.

Definition 2. A continuous linear time-invariant system ẋ(t) = A x(t) + B u(t), y(t) = C x(t) + D u(t), x(0) = x o ,
(2) is said to be positive if for all positive initial states x o i (x o i 0) and all input signals u such that ∀t ∈ R + , u(t) 0, the state vector x(t) and the output vector y(t) are both positive (x(t) 0 and y(t) 0) for all t ∈ R + . Theorem 1. [START_REF] Farina | Positive Linear Systems: Theory and Applications[END_REF]) The linear dynamical system (2) is positive iff A is a Metzler-matrix, B 0, C 0 and D 0.

Estimate of the reachable state of uncertain dynamical systems

Several methods have been developed in the literature to estimate the reachable state space of uncertain systems. In this paper, we have chosen to focus on the one developed in [START_REF] Kieffer | Guaranteed nonlinear state estimation for continuous-time dynamical models from discrete-time measurements[END_REF] and [START_REF] Meslem | Atteignabilité hybride des systèmes dynamiques continus par analyse par intervalles : Application à l'estimation ensembliste[END_REF]. It is based on the theory of differential inequalities and allows to bound the states of an uncertain system by those of two coupled systems without uncertainties.

Consider an uncertain dynamical system described by:

ż = g(z, u, θ), z(0) = z o , (3)
where z(t) ∈ D ⊆ R n is the state vector, u an input signal and θ ∈ [θ ; θ] ⊂ R p a parameter vector of the system ; the vectors θ and θ are the bounds of the parameter vector θ.

The input signal is assumed to belong to the set

V(u, u) = {u : u(t) ∈ [u(t) ; u(t)], ∀t ∈ [0 ; T ]}
, where the signals u and u are both known and satisfy the relation u(t) u(t), ∀t ∈ [0 ; T ]. In the case where l = 1, the Figure 1 below illustrates an example of an input signal u that belongs to the set V(u, u) and an example of an input signal u which does not belong to V(u, u). [START_REF] Kieffer | Guaranteed nonlinear state estimation for continuous-time dynamical models from discrete-time measurements[END_REF] and [START_REF] Meslem | Atteignabilité hybride des systèmes dynamiques continus par analyse par intervalles : Application à l'estimation ensembliste[END_REF], one can determine two vector fields g := g(z, z, u, u, θ, θ) and g := g(z, z, u, u, θ, θ) for which the following two coupled systems ż = g(z, z, u, u, θ, θ) and ż = g(z, z, u, u, θ, θ)

0 T u u u ∈ V(u, u) u ∈ V(u, u) Fig. 1. An illustration of the set V(u, u) in the case where l = 1 Let F ⊆ v ∈ R l : ∃t ∈ [0 ; T ], v ∈ [u(
are such that if z(0) z(0) then, for all t ∈ [0 ; T ] and all triplets (z o , u, θ) ∈ [z(0) ; z(0)] × V(u, u) × [θ ; θ], z(t) z(t) z(t).
The components g i , and g i , i = 1, 2, . . . , n of the vectors fields g and g are constructed as follows:

• g i (z, z, u, u, θ, θ) = g i ξ i , v i , ϑ i , where the function ξ i (resp. v i ) defined on R with values in R n (resp. R l )
and the vector ϑ i ∈ R p are computed as follows:

ξ i j =                  z i if i = j, z j if j = i and ∂g i ∂z j ≥ 0, z j if j = i and ∂g i ∂z j < 0, v i j =          u j if ∂g i ∂u j ≥ 0, u j if ∂g i ∂u j < 0, ϑ i j =          θ j if ∂g i ∂θ j ≥ 0, θ j if ∂g i ∂θ j < 0. • g i (z, z, u, u, θ, θ) = g i ξ i , v i , ϑ i where the function ξ i (resp. v i ) defined on R with values in R n (resp. R l )
and the vector ϑ i ∈ R p are computed as follows:

ξ i j =                  z i if i = j, z j if j = i and ∂g i ∂z j ≥ 0, z j if j = i and ∂g i ∂z j < 0, v i j =          u j if ∂g i ∂u j ≥ 0, u j if ∂g i ∂u j < 0, ϑ i j =          θ j if ∂g i ∂θ j ≥ 0, θ j if ∂g i ∂θ j < 0.

MAIN RESULTS

3.1 An estimate of the reachable output space of the uncertain system S k in (1)

The goal of this section is to give an estimate of the reachable output space of systems S 1 and S 2 . As recalled in the previous section, under some assumptions, the reachable state space of a given uncertain system can be estimated by the state vectors of two specific systems. Some assumptions are made first in this part. Under these assumptions, we show that for each system S k , there exist two LTI systems Σ k and Σ k with outputs y k and y k that are respectively a lower and an upper bound of y k .

Let f k , k = 1, 2 be the two vectors fields associated with the systems S k defined as follows:

f k (x k , u, θ k ) := A k + p k r=1 θ k,r A k,r x k + B k u.
The dependence in t of functions x and u is omitted in the expression of f k for simplicity. Jacobian matrices

∂f k ∂x k = ∂f k i ∂x k,j and ∂f k ∂u = ∂f k i ∂u j
are given by:

∂f k ∂x k = A k + p k r=1 θ k,r A k,r and ∂f k ∂u = B k .
Clearly, signs of

∂f k i ∂x k,j and ∂f k i ∂u j do not change on any subset of R n × R l × R p k .
Therefore, to estimate the reachable state space of systems S k , k = 1, 2 using the method in [START_REF] Kieffer | Guaranteed nonlinear state estimation for continuous-time dynamical models from discrete-time measurements[END_REF] and [START_REF] Meslem | Atteignabilité hybride des systèmes dynamiques continus par analyse par intervalles : Application à l'estimation ensembliste[END_REF], only the signs of the functions ∂f k i ∂θ k,r , i = 1, 2, . . . , n, r = 1, 2, . . . , p k need further consideration. As seen in Proposition 1 below, under Assumptions 1, 2 and 3 on the system S k , these signs do not change on the domain

E k × R l × R p k . Assumption 1.
k ∈ [θ k ; θ k ].
Assumption 2 means that the matrix A(θ k ) is a Metzlermatrix for all parameter vector θ k ∈ [θ k ; θ k ] and that the matrices B k and C k are non-negative. Under Assumptions 1 and 2 on the system S k , without loss of generality, we take E k ⊆ R n + . Indeed, according to the definition of a positive system recalled in Definition 2 and under Assumptions 1 and 2 on the system S k , the vector

x k (t) := x k (t, x o k , u, θ k ) is non-negative for all triplets (x o k , u, θ k ) ∈ X o k × U × [θ k ; θ k ]. Assumption 3. Each row A i k,r , i ∈ {1, 2, . . . , n} of the matrix A k,r = A 1 k,r T A 2 k,r T • • • A n k,r T T is such that A i k,r 0 or A i k,r 0.
Proposition 1. Suppose the systems S k described by (1) satisfy the Assumptions 1, 2 and 3. For all i ∈ {1, 2, . . . , n} and all r ∈ {1, 2, . . . , p k }, the sign of the function

∂f k i ∂θ k,r does not change on the domain E k × R l × R p k . Proof. As E k ⊆ R n + , if A i k,r 0 then, ∂f k i ∂θ k,r (ξ, u, θ) = A i k,r T , ξ ≥ 0, ∀ξ ∈ E k ,
where •, • denotes the usual scalar product on R n . A similar reasoning applied to the case A i k,r 0 implies that

∀ξ ∈ E k , ∂f k i ∂θ k,r (ξ, u, θ k ) ≤ 0.
Lemma 1. Consider the uncertain LTI systems S k described by (1) as follow:

S k ẋk (t) = A k (θ k ) x k (t) + B k u(t), y k (t) = C k x k (t), x k (0) = x o k ∈ X o k . with A k (θ k ) = A k + θ k,1 A k,1 + θ k,2 A k,2 + • • • + θ k,p k A k,p k ,
and suppose the systems S k satisfy the Assumptions 1, 2 and 3. Then, there exist two LTI systems

Σ k    ẋk (t) = A k (θ k , θ k ) x k (t) + B k u(t), y k (t) = C k x k (t), x k (0) = x o k , and
Σ k    ẋk (t) = A k (θ k , θ k ) x k (t) + B k u(t), y k (t) = C k x k (t), x k (0) = x o k , such that the following relation y k (t) y k (t, x o k , u, θ k ) y k (t) holds for all x o k ∈ X o k , all u ∈ U, all θ k ∈ [θ k ; θ k ] and all t ∈ [0 ; T ].
Proof. One can determine with the rules recalled in Section 3 two systems

ẋk (t) = A k (θ k , θ k ) x k (t) + B k u(t), x k (0) = x o k , and
ẋk (t) = A k (θ k , θ k ) x k (t) + B k u(t), x k (0) = x o k , such that the relation x k (t) x k (t) x k (t) holds for all triplets (x o k , u, θ k ) ∈ X o k × U × [θ k ; θ k ] and all t ∈ [0 ; T ]. Therefore, as C k is a non-negative matrix then, C k x k (t) C k x k (t) C k x k (t).
Note that the output vector y k (t) (resp.

y k (t)) of system Σ k (resp. Σ k ) is a function of variables x o k (resp. x o k ), u (resp. u), θ k and θ k i.e y k (t) := y k (t, x o k , u, θ k , θ k ) (resp. y k (t) := y k (t, x o k , u, θ k , θ k )).
The output vectors y k (t) and y k (t) are respectively a lower and an upper bound of y k (t).

Sufficient condition for (X o , U, Θ)-distinguishability

Theorem 2 below gives a sufficient condition for (X o , U, Θ)distinguishability.

Theorem 2. Consider the uncertain LTI systems S k descriebed by (1). Suppose S k satisfy the Assumptions 1, 2 and 3, and consider the output signals

y k (•, x o k , u, θ k , θ k ) and y k (•, x o k , u, θ k , θ k ) of systems Σ k and Σ k obtained in Lemma 1. If there exists an interval I ⊆ [0 ; T ] such that for all t ∈ I, y 2 (t, , x o 2 , u, θ 2 , θ 2 ) ≺ y 1 (t, , x o 1 , u, θ 1 , θ 1 ), (4) or y 1 (t, x o 1 , u, θ 1 , θ 1 ) ≺ y 2 (t, x o 2 , u, θ 2 , θ 2 ), (5)
then, systems S 1 and S 2 are (X o , U, Θ)-distinguishable on [0 ; T ].

Proof. For all (θ 1 , θ 2 ) ∈ Θ and all (x o 1 , x o 2 ) ∈ X,    y 2 (t, x o 2 , u, θ 2 ) y 2 (t, x o 2 , u, θ 2 , θ 2 ), y 1 (t, x o 1 , u, θ 1 , θ 1 ) y 1 (t, x o 1 , u, θ 1 ), ∀t ∈ [0 ; T ] Therefore, if ∀t ∈ I relation (4) holds then y 2 (t, x o 2 , u, θ 2 ) ≺ y 1 (t, x o 1 , u, θ 1 ), ∀t ∈ I. Hence, y 1 (•, x o 1 , u, θ 1 ) = y 2 (•, x o 2 , u, θ 2 ).
If relation (5) holds on I, then the conclusion is the same as previously. It suffices to remark in this case that:

   y 1 (t, x o 1 , u, θ 1 ) y 1 (t, x o 1 , u, θ 1 , θ 1 ), y 2 (t, x o 2 , u, θ 2 , θ 2 ) y 2 (t, x o 2 , u, θ 2 ), ∀t ∈ [0 ; T ].
The sufficient condition for (X o , U, Θ)-distinguishability of Theorem 2 concerns with positive LTI systems, which constitute an important class of LTI systems [START_REF] Farina | Positive Linear Systems: Theory and Applications[END_REF]) and which includes in particular many physical systems (networks of reservoirs, heat exchangers and distillation columns,...). If this condition is not satisfied, one cannot draw any conclusion on the distinguishability/indistinguishability between S 1 and S 2 . The two tanks have both the same section σ = 1. Assume that the flows d i , i = 1, 2 of the fluid are such that

An academic example

d i = a i h i where a i = a o i +α i , α i ∈ [α i ; α i ] and a o i -α i ≥ 0. Choosing x 1 (t) = h 1 (t) h 2 (t) and y 1 (t) = h 2 (t).
Therefore, a state space representation of S 1 is

S 1 :        ẋ1 (t) = (A 1 + α 1 A 1,1 + α 2 A 1,2 ) x 1 (t) + 1 0 u(t), y 1 (t) = [0 1] x 1 (t), x 1 (0) = x o 1 ∈ [x o 1 ; x o 1 ]
, where

A 1 = -a o 1 0 a o 1 -a o 2 , A 1,1 = -1 0 1 0 and A 1,2 = 0 0 0 -1 . For all (α 1 , α 2 ) ∈ [α 1 ; α 1 ] × [α 2 ; α 2 ], the matrix A 1 (α) = A 1 + α 1 A 1,1 + α 2 A 1,2 with α = α 1 α 2 is a
Metzler-matrix and matrices A 1,i , i = 1, 2 satisfied the conditions of Assumption 3. A simple application of the computation rules of Σ 1 and Σ 1 gives:

Σ 1 :        ẋ1 (t) = -a o 1 -α 1 0 a o 1 + α 1 -a o 2 -α 2 x 1 (t) + 1 0 u(t), y 1 (t) = [0 1] x 1 (t), x 1 (0) = x o 1 , and
Σ 1 :        ẋ1 (t) = -a o 1 -α 1 0 a o 1 + α 1 -a o 2 -α 2 x 1 (t) + 1 0 u(t), y 1 (t) = [0 1] x(t), x(0) = x o 1 .
In addition to system S 1 , consider the two tanks system S 2 of Figure 3 below. Assume that

• the two tanks have both the same section σ = 1.

• the flows di , i = 1, 2 and d12 are such that di = b i hi , i = i, 2 and d12 = b 12 h2 -h1 .

• Therefore, a state space representation of S 2 is

S 2 :        ẋ2 (t) = (A 2 + β A 2,1 ) x 2 (t) + 0 1 u(t), y 2 (t) = [0 1] x 2 (t), x 2 (0) = x o 2 ∈ [x o 2 ; x o 2 ]
, where

A 2 = -b o 2 -b 12 b 12 b 12 -b 12 -b 1 and A 2,1 = -1 0 0 0 .
The systems Σ 2 and Σ 2 which bound S 2 are defined as follows:

Σ 2 :        ẋ2 (t) = A 2 (β) x 2 (t) + 0 1 u(t), y 2 (t) = [0 1] x 2 (t), x 2 (0) = x o 2 and Σ 2 :        ẋ2 (t) = A 2 (β) x 2 (t) + 0 1 u(t), y 2 (t) = [0 1] x 2 (t), x 2 (0) = x o 2 , where A 2 (β) = -b o 2 -b 12 -β b 12 b 12 -b 12 -b 1 and A 2 (β) = -b o 2 -b 12 -β b 12 b 12 -b 12 -b 1
The numerical results of the Figure 5 are obtained in the case where

• a o 1 = 0.1 = a o 2 , α 1 = α 2 = -0.005, α 1 = α 2 = 0.005. • b 1 = b o 2 = 0.2, b 21 = 0.1, β = -0.005, β = 0.005. • x o 1 = 0 0 = x o 2 , x o 1 = 0.02 0.02 = x o 2 , T = 500.
• the space U = V(u, u) of admissible inputs of the systems S 1 and S 2 is defined by the piecewise constant signals u and u of the Figure 4 below. • U is defined by the Figure 4 and Θ = [-0.005 ; 0.005] 3 .

CONCLUSIONS AND OPEN PROBLEMS

Distinguishability is a property that allows one to determine the active mode of a switched system by only observing its input-output data. In the literature, most of the works dealing with the subject only deal with LTI systems without uncertainties. This paper proposes to deal with linear systems with affine uncertainties in the state space model. In this new context, an adapted definition of distinguishability is given. Two systems can be distinguished when their reachable output spaces are disjoint.

When the considered systems are positive, extension of the result in [START_REF] Kieffer | Guaranteed nonlinear state estimation for continuous-time dynamical models from discrete-time measurements[END_REF] and [START_REF] Meslem | Atteignabilité hybride des systèmes dynamiques continus par analyse par intervalles : Application à l'estimation ensembliste[END_REF] allows to bound all their admissible outputs by the outputs of two specific LTI systems. An estimate of the reachable output space of these systems in response to an input belonging to an uncertain range can then be determined. A sufficient condition for distinguishability is derived from this estimation.

As the sufficient condition that we have obtained strongly depends on the estimation of the reachable output spaces, a further work will aim to find a method to obtain a better estimate of the reachable output space of the considered systems. On the other hand, we have shown in [START_REF] Motchon | Operating modes distinguishability condition in switching systems[END_REF] that, when uncertainties are ignored in the systems' models, the conditions of distinguishability can be deduced from some of intrinsic properties of the models.

In the same way, for the class of linear systems with affine uncertainties, it would be interesting to determine the intrinsic properties which can ensure the distinguishability of two systems.

Figure 2 Fig. 2 .

 22 Figure 2 below illustrates a system S 1 of two tanks. P

Fig. 4 .Fig. 5 .

 45 Fig. 4. Space U of admissible inputs of systems S 1 and S 2

  b 2 = b o 2 + β where β ∈ [β ; β].

		u
	P	h1	d12	h2
		d1		d2
	Fig. 3. System S 2 of two tanks
	Choosing
	x 2 (t) = h2 (t) h1 (t)	and y 2 (t) = h1 (t).