Guillaume Bouleux
email: guillaume.bouleux@univ-st-etienne.fr

Oblique Projection pre-Processing and TLS application for diagnosing rotor bar defects by improving power spectrum estimation

Keywords: Induction Machine, fault diagnosis, Motor Current Signal Analysis, frequency estimation, rotor bars

Diagnosing defects on rotating machines can be reached by several angles. When dealing with asynchronous motor drive, such physical elements rotate that a natural angle for treating the healthiness of the motor can be obtained by the use of spectral analysis tools. It is now stated that electrical or mechanical defects, which appear periodically as well, can be retrieved by analyzing the amplitude of particular frequencies inside an estimated power spectrum. When dealing with broken rotor bars detection it is essential to accurately localize the frequencies related to the slip inside the power spectrum. The diagnosis is thereafter made by indicators given with respect to their power. For actual low level of load operations, the supply frequency generally masks the frequencies which could be exploited for indicators. Therefore, we propose to cancel, as well as possible, the contribution of this supply frequency to develop the useful and closely components. The resolution should be very thin for the components to be estimated. In consequence, we use a prior-knowledge subspace-based frequency estimator, already developed in the literature, we complete with an Oblique Projection coupled with a Total Least Squares solution for estimating the power of the resulting estimated frequencies. Finally, we show by means of a real application how it contributes to improve the power spectrum estimation when compared to the FFT or periodogram-based analysis and how the aforementioned power spectrum makes the diagnosis indicator of rotor bars efficient.

Introduction

It is a fact, the diagnosis of electrical or mechanical defects of an induction motor has concerned and focused the attention since the last two decades [START_REF] Benbouzid | A Review of Induction Motors Signature Analysis as a Medium for Faults Detection[END_REF][START_REF] Jung | Online diagnosis of induction motors using MCSA[END_REF][START_REF] Menacer | Experimental Detextion of Rotor Failures of Induction Machines by Stator Current Spectrum Analysis in Function of the Broken rotor Bars Position and the Load[END_REF][START_REF] Peng | Application of the wavelet transform in machine condition monitoring and fault diagnostics: a review with bibliography[END_REF][START_REF] Jardine | A review on machinery diagnostics and prognostics implementing condition-based maintenance[END_REF][START_REF] Jayaswal | Review Article Machine Fault Signature Analysis[END_REF][START_REF] Mehrjou | Rotor fault condition monitoring techniques for squirrel-cage induction machine-a review[END_REF][START_REF] García-Escudero | Robust condition monitoring for early detection of broken rotor bars in induction motors[END_REF][START_REF] Tran | An intelligent condition-based maintenance platform for rotating machinery[END_REF][START_REF] Didier | A new approach to detect broken rotor bars in induction machines by current spectrum analysis[END_REF][START_REF] Zhang | A frequency domain approach to fault detection in sampled-data systems[END_REF][START_REF] Dalpiaz | Condition monitoring and diagnostics in automatic machines: Comparison of vibration analysis techniques[END_REF]. Based on the electrical modeling of the machine [START_REF] Abed | Modelling Induction Motors for Diagnostic Purposes[END_REF][START_REF] Casimir | Comparison of modelling methods and of diagnostic of asynchronous motor in case of defects[END_REF] and on the choice of the sensor-type, several approaches have been proposed [START_REF] Cristalli | Mechanical fault detection of electric motors by laser vibrometer and accelerometer measurements[END_REF][START_REF] Siddique | A review of stator fault monitoring techniques of induction motors[END_REF][START_REF] Worden | Natural computing for mechanical systems research: A tutorial overview[END_REF][START_REF] Widodo | Support vector machine in machine condition monitoring and fault diagnosis[END_REF][START_REF] Lei | A review on empirical mode decomposition in fault diagnosis of rotating machinery[END_REF]. Usually, but not only, the major manner of acquiring a significant signal from the motor is made by acoustic [START_REF] Hamzaoui | Vibro-acoustic analysis and identification of defects in rotating machinery, part i: Theoretical model[END_REF], vibratory [START_REF] Randall | Vibration-based Condition Monitoring: Industrial, Aerospace and Automotive Applications[END_REF][START_REF] Didier | Fault detection of broken rotor bars in induction motor using a global fault index[END_REF] or current acquisition [START_REF] Eltabach | Broken rotor bars detection by a new non-invasive diagnostic procedure[END_REF][START_REF] Eltabach | Quantitative analysis of noninvasive diagnostic procedures for induction motor drives[END_REF]. Each has its advantages and of course its drawbacks; it will depend on the default characterization and eventually the motor drive particularity [START_REF] Randall | Vibration-based Condition Monitoring: Industrial, Aerospace and Automotive Applications[END_REF][START_REF] Bellini | Advances in Diagnostic Techniques for Induction Machines[END_REF][START_REF] Silva | Induction motor fault diagnostic and monitoring methods[END_REF][START_REF] Elbouchikhi | Induction machine fault detection enhancement using a stator current high resolution spectrum[END_REF]. However, we point out that all these technics measure periodic information directly came from the motor. The inherent rotation of bearings, magnetic fields or currents inside the motor, is relevant information that a vibratory, an acoustic or obviously a current signal carry on. Accepting the time periodicity allows us to relate the different faults to a frequency-based examination of the signal acquired. Depending on the severity of the default, the power of the related frequencies will clearly evolve. The diagnosis can in consequence be performed by analyzing an indicator based on the power ratio between the supply frequency and its sidebands components [START_REF] Bellini | Advances in Diagnostic Techniques for Induction Machines[END_REF][START_REF] Elbouchikhi | Induction machine fault detection enhancement using a stator current high resolution spectrum[END_REF]; the more this ratio increases the more faulty the motor is.

The first works on this topic were focused on Power Spectral Density (PSD) or related and closely works essentially based on periodogram-like treatment. The non-parametric estimation procedure gives then estimation of both, frequencies and amplitudes needed for obtaining an efficient diagnostic indicator. These kinds of approaches were certainly first confidence approach giving reliable results. However, for broken rotor bars detection, an excellent accuracy location of the sidebands components is essential [START_REF] Mehrjou | Rotor fault condition monitoring techniques for squirrel-cage induction machine-a review[END_REF][START_REF] King | Current signature analysis of induction machine rotor faults using the fast orthogonal search algorithm[END_REF]. Then, overcoming the resolution insufficiency of periodogram-based approaches which we remind, is limited by the Fourier resolution, i.e. the ratio between the sampling frequency and the numbers of samples; high resolution methods have emerged, e.g. the Fast Orthogonal Search (FOS) algorithm [START_REF] Mcgaughey | Speed sensorless estimation of ac induction motors using the fast orthogonal search algorithm[END_REF][START_REF] Chon | Accurate identification of periodic oscillations buried in white or colored noise using fast orthogonal search[END_REF], some MUSIC-like method [START_REF] Elbouchikhi | Induction machine fault detection enhancement using a stator current high resolution spectrum[END_REF][START_REF] Bracale | A High Resolution Method for On Line Diagnosis of Induction Motors Faults[END_REF][START_REF] Kia | A High-resolution frequency estimation method for three-phase induction machine fault detection[END_REF][START_REF] Bouleux | A subspace-based rejection method for detecting bearing fault in asynchronous motor[END_REF] or ESRIT-like improvement [START_REF] Bouleux | An improper random vector approach for esprit and unitary esprit frequency estimation[END_REF] and Markov-based methods for frequency estimation [START_REF] Wirfalt | Subspace-based frequency estimation utilizing prior information[END_REF]. The high resolution property of these methods is obtained thanks to an algebraic parameterization of the acquisition signal which permits the use of mathematical concepts of vector subspaces. Two estimation steps are thus required to fulfill the spectral analysis. Let us summarize the steps by the following: (i) the frequencies are obtained when the subspaces, with respect to those parameters, fit the true ones as closely as possible. This notion gives rise to the more general concept of subspace fitting [START_REF] Kristensson | Further results and insights on subspace based sinusoidal frequency estimation[END_REF] which group together a large number of subspace-based solutions, including Maximum Likelihood (ML) frequency estimation;

(ii) a classical ordinary Least Squares (LS) procedure is utilized for obtaining the power of each estimated frequencies.

Owing to its non-invasive or non-intrusive capability [START_REF] Eltabach | Quantitative analysis of noninvasive diagnostic procedures for induction motor drives[END_REF], stator current acquisition is increasingly put in place. For squirrel-cage rotor induction motors, electrical defects such as rotor bar defects could be seen inside the stator current signal by monitoring the frequencies around the supply frequency [START_REF] Benbouzid | A Review of Induction Motors Signature Analysis as a Medium for Faults Detection[END_REF][START_REF] Mehrjou | Rotor fault condition monitoring techniques for squirrel-cage induction machine-a review[END_REF][START_REF] King | Current signature analysis of induction machine rotor faults using the fast orthogonal search algorithm[END_REF]. In such an environment, the driven frequency is undisputedly the most powerful frequency inside the stator current spectrum and it may mask the meaningful sidebands frequencies. To tackle this problem, recent work incorporated the knowledge of the driven frequency into subspace methods for improving the estimation of the remaining frequencies and getting rid of the known one [START_REF] Bouleux | A subspace-based rejection method for detecting bearing fault in asynchronous motor[END_REF][START_REF] Wirfalt | Subspace-based frequency estimation utilizing prior information[END_REF][START_REF] Kristensson | Further results and insights on subspace based sinusoidal frequency estimation[END_REF]. The question of improving the power estimation of the frequencies that have not been eliminated, and only those ones was an open question to which we respond in the present paper.

The paper is organized as follows. We first state the model of observation in section 2, we then recall in the third section the subspace estimation procedure which takes into account the knowledge of some frequencies; section 4 presents the Oblique projection coupled to the Total Least Squares (TLS) solution for improving power spectrum estimation and finally we analyze the method through an application in section 5.

Model Specification

A classical Hilbert transform is first assumed for suppressing negative frequencies and avoiding aliasing effects. We assume therefore that the signal acquired is the superposition of d complex sinusoids corrupted by an acquisition noise. We suppose the temporal acquisition of N samples of the form

z(t) = d k=1 α k e j(ω k +φ k ) + n(t), t = 1, . . . , N (1)
where α k is a deterministic and constant amplitude, ω k is the deterministic pulsation (we will mix up the pulsation and the frequency here and after), φ k is the initial phase at t = 0 and assumed to be a random variable taking its values uniformly in the range [0; 2π(. We also assume that n(t) is a random circular and zero mean Gaussian variable independent of {ω k } d k=1 with covariance given by

E[n(t 1 )n ⋆ (t 2 )] = σ 2 δ(t 1 , t 2 ), (2)
where E[.] stands for the mathematical expectation, the superscript . ⋆ corresponds to the complex conjugate and where δ(t 1 , t 2 ) = 1 for t 1 = t 2 and is equal to zero for any other values. In order to use a subspace approach, whereas only one sensor is used, we define a data vector y(t) viewed as a windowed partition of the whole data set acquired. Then, y(t) = [z(t) . . . z(t + m -1)] T , with m the window length and . T the transpose operator. Henceforth, bold is used to specify the vectors and bold capitals is used for matrix distinction. Consequently, we obtain the so-called matrix representation

y(t) = ASw(t) + n(t) (3)
in which

A =      1 . . . 1 . . . . . . e jω1(m-1) . . . e jω d (m-1)      , S =      α 1 e j(φ1) 0 . . . 0 α d e j(φ d )      , (4)
and where w(t) = e jω1t . . . e jω d t . The complex vector n(t) is a random circular and zero mean Gaussian vector with covariance E[n(t)n H (t)] = σ 2 I where . H means conjugate transpose and where I is the identity matrix. The algebraic representation of the linear application between the input complex sinusoids and the (Nm + 1) virtual output realization of the stochastic process is reached by staking the (Nm + 1)

windowed data vectors such as

Y = y(t) . . . y(t + N -m + 1) . (5)
For an initial time t = 0, we obtain the matrix Y = [y(0) . . . y(Nm + 1)] which is a Hankel matrix whose the basic form is

Y =         y(0) y (
        = ASW + N , (6)
with W = [w(0) . . . w(Nm + 1)] and N = [n(0) . . . n(Nm + 1)]. Then, the output matrix can be built by staking in column-wise the (Nm + 1) data windowed samples of the original and available data set.

For ease of reading and notation, we have not explicitly written the ω dependence of A. It is however necessary to keep in mind that the frequencies will be estimated thanks to the subspace spanned by the columns of A and that the amplitudes and a posteriori the power of the sinusoids will be estimated thanks to the S matrix. The challenge is thus to estimate both the frequencies {ω k } d k=1 and the amplitudes {α k } d k=1 with only the observation matrix Y .

In modeling the asynchronous motor drive electrically, we can easily point out the relations between a default on one or several rotor bars and the electromotive force of the machine (cf. references inside [START_REF] Bellini | Advances in Diagnostic Techniques for Induction Machines[END_REF]). Indeed, the inductive current admits an analytic expression in which clearly the defaults are visible through the apparition of harmonics around the on line frequency or the driven frequency. Basically, we have the well known equality

f sb = (1 + 2kg)f r , k ∈ Z *
where g is the slip coefficient, f r the driven frequency which could be 50Hz or 60 Hz and f sb the frequency of sidebands components. For monitoring as well as possible and diagnose with a high confidence level, it is sufficient to at least inspect the sidebands components; but this could be made difficult if the frequencies are close to the online one, especially for low level of load operations. In addition, the power of, for example, the 50 Hz, is by far the most important one. These conditions pushed up for treating the problem with a new viewpoint. The interest is no more the estimation of the whole frequency set {ω k } d k=1 but only a subset where the 50 Hz, and possibly its harmonics, would be somehow mathematically filtered but useful for improving the estimation of the other frequencies. The challenge is in consequence to utilize the 50 Hz (and its harmonics if any) to improve the estimation of a subset {ω k } du k=1 , with d u the number of unknown frequencies, and to estimate thereafter their associated power. The two following sections are dedicated to these tasks.

Estimating the frequencies with support of 50 Hz

We consider that there already exist solutions to this problem, this section will then draw the principal steps of the method; to have a comprehensive idea, we advise the reader to refer themselves to the works [START_REF] Wirfalt | Subspace-based frequency estimation utilizing prior information[END_REF][START_REF] Bouleux | An Optimal Prior knowledge-based DOA Estimation Method[END_REF] and the references therein.

The method principle

Let start by introducing the following polynomial

Q(z) b 0 z d + b 1 z d-1 + . . . + b d = b 0 Π d k=1 (z -e jω k ) (7)
for which the argument of each zero is the frequency of interest. It is easy to incorporate the knowledge of the 50 Hz and its harmonics inside this polynomial. It suffices for this to factorize Q(z) by the known zeros associated with the known frequencies. We obtain consequently the following expression

Q(z) = P u (z)P k (z) (8)
where the polynomial P k (z) is composed by the known d k zeros/frequencies and the polynomial P u (z) by the d u unknown zeros/frequencies in the following manner

P u (z) b0 z du + b1 z du-1 + . . . + bdu = b0 Π du k=1 (z -e jω k ) (9)
P k (z) c 0 z d k + c 1 z d k -1 + . . . + c d k = c 0 Π d k k=1 (z -e jω k ). (10)
The problem of frequency estimation is then reported to the estimation of the unknown polynomial coefficients { bk } du k=1 .

How to reach the frequencies

Owing to the definition of Q(z), we have the crucial algebraic property

B H A = 0 (11)
in which the matrix B H of size (md) admits the structure

B H =         b d . . . b 0 0 . . . 0 0 b d . . . b 0 0 . . . . . . 0 . . . . . . . . . 0 0 . . . 0 b d . . . b 0         (12)
This relation basically means that the vector b H = [b d . . . b 0 ] belongs to the kernel of the A matrix; if we assume that the columns of A span the signal subspace then we will say that b belongs to the noise subspace in reference to the statistical independence and representation and that of noise. The set {b k } d k=1 could be consequently estimated when the orthogonality between b and the columns of A is maximal. It is then necessary to have access to the signal subspace or an estimate of it.

Let start with the sample covariance matrix defined by

R = 1 N -m + 1 N -m+1 t=1 y(t)y H (t) (13)
and perform next an eigenvalue decomposition of the sample covariance to obtain

R = Ês Λs ÊH s + Ên Λn ÊH n ( 14
)
where Define next êsk as the k-th vector of Ês . A Maximum Likelihood (ML) estimation procedure is performed from the random sequence {B H ês1 , . . . , B H êsd } and after algebraic manipulations, b is easily obtained by minimizing the quadratic form b H Qb where Q is explicitly given in reference [START_REF] Bouleux | An Optimal Prior knowledge-based DOA Estimation Method[END_REF] where the reader will find all materials needed for the criterion derivation.

We now use the convolution property

b = C b (15)
with

C T =         c 0 . . . c d k 0 . . . 0 0 c 0 . . . c d k 0 . . . . . . 0 . . . . . . . . . 0 0 . . . 0 c 0 . . . c d k         (16)
and where C is of size (d + 1) × (d u + 1), to give below, the minimization criterion from which the { bk } du k=1 will be estimated

min b s.t. b 2 bH C T QC ⋆ b ( 18
)
with . ⋆ the conjugate operator. Different hints for solving (18) can be found in [START_REF] Wirfalt | Subspace-based frequency estimation utilizing prior information[END_REF][START_REF] Kristensson | Further results and insights on subspace based sinusoidal frequency estimation[END_REF][START_REF] Bouleux | An Optimal Prior knowledge-based DOA Estimation Method[END_REF]. Once b is obtained, the frequencies are deduced from the roots of P u (z).

We have now to treat the problem of the amplitudes/power for fulfilling a good diagnosis.

The TLS-based amplitude estimation

When the frequencies have been estimated the acquisition model ( 6) is modified to explicitly point out the dependence on the known and unknown frequency contribution and it yields

Y = [A(ω k ) A(ω u )] Â(ω)   S k S u   S W k W u W +N ( 19
)
where the Vandermonde matrix Â(ω) has been portioned between the sub-matrix A(ω u ) of estimated unknown frequencies and the sub-matrix A(ω k ) of known frequencies. Based on [START_REF] Lei | A review on empirical mode decomposition in fault diagnosis of rotating machinery[END_REF] we aim at estimating only the unknown amplitudes. For this to be achieved, we have to first cancel off the contribution of known quantities such as the frequencies and the amplitudes from ( 19) and then only estimate by a Total Least Squares (TLS) approach the remaining amplitudes. This is summed up by 1. Obtaining a model which could be similar to

Y r = A(ω u )S u W u + N r (20)
2. Estimating the amplitudes Ŝu from (20).

This two stage procedure is developed in the next sections.

Preprocessing the data observation

Reducing the number of linear combinations inside ( 19) is a task which has been of some concern in several works [START_REF] Chen | Parameter Estimation with Prior Knowledge of Known Signal Poles for the Quantification of NMR Spectroscopy Data in the Time Domain[END_REF]. A possible way to tackle this problem is to basically null the known column-subspace A(ω k ) while leaving the unknown column-subspace A(ω u ) unaltered. The work of Behrens and Scharf [START_REF] Behrens | Signal Processing Applications of Oblique Projection Operators[END_REF],

we briefly recall, is a lead for solving this problem.

One can start by examining the orthogonal projection onto the entire signal subspace when it explicitly depends on the known and unknown subspaces. We have then

P [A(ω k ) A( ωu)] = [A(ω k ) A(ω u )]     A H (ω k ) A H (ω u )   A(ω k ) A(ω u )   -1   A H (ω k ) A H (ω u )   = [A(ω k ) A(ω u )]     A H (ω k )A(ω k ) A H (ω k )A(ω u ) A H (ω u )A(ω k ) A H (ω u )A(ω u )     -1   A H (ω k ) A H (ω u )   (21)
using the inversion matrix lemma, we straightforwardly obtain that

    A H (ω k )A(ω k ) A H (ω k )A(ω u ) A H (ω u )A(ω k ) A H (ω u )A(ω u )     -1   A H (ω k ) A H (ω u )   =   (A H k P ⊥ Âu A k ) -1 A H k P ⊥ Âu ( ÂH u P ⊥ A k Âu ) -1 ÂH u P ⊥ A k   (22)
where for ease of notation A(ω k ) and A(ω u ) have been replaced by A k and Âu respectively, P ⊥ A k is the orthogonal projection onto the complementary subspace of A k (i.e. the null space N (A k )) basically defined by

P ⊥ A k = I -A k A H k A k -1
and similarly P ⊥ Au = Âu ÂH u Âu -1 ÂH u is the orthogonal projector onto N ( Âu ). It clearly follows from [START_REF] Didier | Fault detection of broken rotor bars in induction motor using a global fault index[END_REF] that the orthogonal projector onto the entire signal subspace decompose itself as

P [A(ω k ) A( ωu)] = P (A k | Âu) + P ( Âu|Ak) . (23)
The projection operators P (A k | Âu) and P ( Âu|Ak) admit the following definitions

P (A k | Âu) = A k (A H k P ⊥ Âu A k ) -1 A H k P ⊥ Âu ( 24
)
P ( Âu|Ak) = Âu ( ÂH u P ⊥ A k Âu ) -1 ÂH u P ⊥ A k (25)
and are seen as oblique projectors. Indeed, P ( Âu|Ak) projects onto the subspace spanned by the columns of Âu , along with the subspace spanned by the columns of A k . The effects of the vectors of A k are therefore annulled while the information contained in the columns of Âu is unaffected. For more explanations and geometrical interpretations we invite the reader to see [START_REF] Behrens | Signal Processing Applications of Oblique Projection Operators[END_REF]. In conclusion, we build the oblique projector P ( Âu|Ak) from the estimated frequencies and the known frequencies (50 Hz and its harmonics for example), and project the data observation matrix Y onto the unknown subspace spanned by the columns of Âu along with the known subspace. This way, we now have an observation model equivalent to [START_REF] Hamzaoui | Vibro-acoustic analysis and identification of defects in rotating machinery, part i: Theoretical model[END_REF]. The amplitudes can in consequence be estimated via general linear estimation procedures.

The TLS procedure for amplitude estimation

In this section, we propose to overcome the difficulty of estimating the amplitudes. Our approach differs from general Least Squares (LS) fit to the data observation in a sense that (i) we use the set of Nm + 1 vectors instead of using directly the N available samples and (ii) we further use a TLS fit.

If we define

y r0 = P ( Âu 0 |A k 0 ) y 0 = Âu0 s u0 + n 0 (26)
with y 0 the vector of N available data observation, Âu0 the Vandermonde matrix Âu extended to N rows,

s u0 = [α 1 e φ1 .
. . α du e jφ du ] T the vector of complex amplitudes and n 0 the N samples additive white Gaussian noise vector; the parameter estimation is generally solved by employing a LS procedure for obtaining ŝu0 = Â♯ u0 y r0 [START_REF] Elbouchikhi | Induction machine fault detection enhancement using a stator current high resolution spectrum[END_REF] in which . ♯ denotes the Moore-Penrose pseudo-inverse. Finally, the power of the frequencies is deduced by taking the 2-norm of the estimated vector ŝu0 . We prefer here to use the stochastic notion of model (20) (or originaly ( 6)) to directly estimate the amplitudes by using the Nm + 1 stochastic realizations. Our model is consequently [START_REF] Hamzaoui | Vibro-acoustic analysis and identification of defects in rotating machinery, part i: Theoretical model[END_REF] and no more [START_REF] Silva | Induction motor fault diagnostic and monitoring methods[END_REF]. In addition, since Âu has been estimated, we can consider, without loss of generality, that the columns of Âu have some erroneous values; if we admit that the observations Y r = P ( Âu|Ak) Y are measured with some uncertainty, which once again is sufficiently realistic, we fall into the multidimensional TLS problem and not the LS one.

The multidimensional TLS problem [START_REF] Huffel | The total least squares problem : computational aspects and analysis[END_REF] minimizes at the same time the errors of both Y r and Âu in establishing the relation

[ Âu Y r ]   S u W u -I   ≈ 0 (28)
where the approx ≈ symbol reflects the error on Y r (the white Gaussian noise N r ). This means that

[ Âu Y r ] = U Σ   V 11 V 12 V 21 V 22   (29)
in which [V T 11 V T 21 ] T is of size (m×d u ) and thus resolves the best d u -rank approximation while the remaining (Nmd u+1 ) right eigenvectors, i.e. [V T 12 V T 22 ] T , resolve the null space of [ Âu Y r ]. The SVD decomposition is not unique and therefore there exists a unitary permutation matrix Ξ such as

  S u W u -I   =   V 12 V 22   Ξ. (30)
From [START_REF] Chon | Accurate identification of periodic oscillations buried in white or colored noise using fast orthogonal search[END_REF] we have first the relation

Ξ = -V -1 22 (31)
from which we deduce the TLS amplitudes parameter estimation

Ŝu Ŵ u = -V 12 V -1 22 . (32)
We mean in fine over the (Nm + 1) realizations of the amplitudes for having the final estimation. The power of each sinusoid is thereafter computed by taking the 2-norm of the corresponding vector.

With the two previous main sections, we have given the tools for estimating the frequencies and their associated power for monitoring and diagnosing rotor bars defects from a noisy stator current measurement.

The following section is dedicated to the analysis of the proposed methods in a real life scenario.

Diagnosing rotor bar defects

This section focuses on the application of the methods proposed on a real asynchronous motor for diagnosing the case of 1 or 2 broken bars on the squirrel-cage rotor by means of a power spectrum analysis.

We first present the test bench and the experimental conditions and next we analyze the different results

given by the application of the method.

Test bench and experimental conditions

An asynchronous motor drive of 4KW with a squirrel-cage rotor was used for the experiments. The AC machine was coupled with a DC motor with 2 pole pairs and manufactured by the Entreprise Algerienne des Industries Electrotechniques Electro-Industries-AZAZGA company. We have measured three rotor stateshealthy, 1 and 2 broken bars. Fig. 1 illustrates the cases of healthy rotor and 2 broken bars conditions. The purely subjective choice of drilling rotor bars at the junction of one bar and the shorted ring was adopted due to the high failure rate in this position. Indeed, it is more likely that the bars will break at the welding between both elements. Finally, the load varies thanks to resistances connected to the DC drive which in turn becomes a DC generator.

We have used hall effect sensors for measuring voltage and current signals and they have been acquired via an acquisition card IOTEQ/DAQ 1005 series with a sampling frequency f s =200KHz, 16 analog inputs and 16 digital inputs.

Improving the Power Spectrum estimation

To be in accordance with the complex data model observation ( 6 Actually, accuracy is reached but resolution is a bit lost with this method.

We then wanted to have an idea whether a refined FFT-algorithm, such as the ZOOMFFT technique (see the references of [START_REF] Mehrjou | Rotor fault condition monitoring techniques for squirrel-cage induction machine-a review[END_REF]) could improve the estimation of the sidebands. The principle is to low-pass filter the 10s of the healthy machine at nominal level of load with a cut of frequency of 100 Hz (Fig. 3-(a)) and next downsample it by a factor of 10. We observe that once again, this procedure is unable to estimate the sidebands components around the driven frequency, Fig. 4 draws that fact.

In short, we have shown that for 10s of acquisition, low-pass filtered and downsampled by a factor of 10, leading to a new sampling frequency equal to 1 KHz; the proposed method which integrates the 50 Hz and for which an oblique pre-processing followed by a TLS solution is employed was clearly better than classical PSD or FFT-based analysis for spectrum estimation. This fundamental observation, developed for the case of a healthy machine at nominal level of load is furthermore reinforced by the results depicted in Fig. 5 and Fig. 6 for which the machine has 1 broken bar at respectively nominal and low level of load. Finally, with no ambiguity, whatever the load and the number of broken bars, Fig. 5 to Fig. 8, the proposed method undisputedly gives the frequencies related to the sidebands. This estimation will next help us when it comes to using a diagnosis indicator.

A diagnosis though improved

In this part, we exploit the frequencies related to the sidebands components for they help to have an indicator of the rotor bars state. We use so the concept of power ratio between the supply frequency and the sidebands components [START_REF] Bellini | Advances in Diagnostic Techniques for Induction Machines[END_REF], with for example [START_REF] Elbouchikhi | Induction machine fault detection enhancement using a stator current high resolution spectrum[END_REF]

I = k∈S s 2 k s 2 0 where S = {k ∈ Z * | f sb = (1 + 2kg)f r is estimated}.
We have used this indicator in conjunction with the high resolution method proposed in [START_REF] Elbouchikhi | Induction machine fault detection enhancement using a stator current high resolution spectrum[END_REF] for the same conditions as for the experience of Fig. 3-(b), 5 and 7.

The results are reported inside table1. If we only inspect the indicator values of table 1, we find no ambiguity between healthy and damaged rotor bars state. Clearly, we have a ratio between healthy and 1 broken bar condition, approaching the value of 200 and a ratio between 1 and 2 broken bars around 6. However, the indicator is based on wrong estimated frequencies. The supply frequency is not accurately estimated by the method for at least broken bars scenario and sidebands components are in total error as well, even if their number could be reasonably coherent. In fact, the method of [START_REF] Elbouchikhi | Induction machine fault detection enhancement using a stator current high resolution spectrum[END_REF] is very sensitive to noisy conditions and the fact the frequencies are estimated by a fixed structure gives a chance to estimate almost everything, leaving a room for artifacts.

In conclusion what is the importance/confidence we can give to an indicator built from wrong information.

The diagnosis is in consequence questioned and the decision finally becomes hard to make.

Based on the estimated power spectrum of the method proposed in this paper (Fig. 5 to Fig. Whatever the load conditions, thanks to a better estimated power spectrum, the sidebands components are better estimated and the diagnosis indicator is more significant . From the results of tables 2 and table 3, we easily see the difference between healthy and failure conditions. The ratio between a healthy and damaged condition is around 30, which leaves no doubt about the decision to make. Similarly, the difference between 1 and 2 broken bars is significant. We have first the number of harmonics growing from 1 to 2 with odd harmonics accurately estimated and next inevitably the increase of the indicator value. The ratio between 1 and 2 broken bars of the latter is for example in between 4 and 7 with the highest difference for the low level of load.

Conclusion

One of the ways to come to a reliable diagnosis on a machine from current acquisition is to inspect an indicator based on the ratio between the supply frequency and its sidebands components. For different actual operations, the load strongly impacts on the proximity of the sidebands to the driven frequency. The challenge is therefore to be able to estimate these sidebands and especially for low level of load. Classical approaches based on periodogram or FFT-based analysis have shown their inadequacy for such a difficult problem. The idea developed in this paper is then to overcome this difficulty by using a high resolution method designed with the support of the supply frequency. The purpose is therefore to suppress this particular and harmful frequency in the set of frequencies of interest. This method has already been introduced and it gives accurate frequency estimation with a sufficient resolution, allowing to precisely discriminate frequencies such as sidebands. The major point is that the driven frequency no longer exists. The next treatment consists in estimating the power associated with the remaining frequencies. In this paper, we propose to first pre-process the data observation by obliquely projecting it onto a subspace which would be associated with the spectrum of interest, i.e. without the supply frequency, along the subspace in which the supply frequency would live. This consequently accounts for the supply frequency in the data observation and then cancels (almost) its influence on the power spectrum. Once this pre-processing of the data observation has been done, we propose a Total Least Squares (TLS) estimation procedure to finally have the power of each frequency without the impact of the very powerful supply frequency. We have tested this pro-cedure in different conditions of load, nominal or low level, for healthy, 1 broken and 2 broken bars against classical FFT analysis, ZOOMFFT method and MUSIC-like criterion especially designed for mechanical defects diagnosing. The results observed for short acquisition duration have clearly shown the limits of these methods. Conversely, the diagnosis given by an indicator based on the power spectrum estimated by the present approach is sure and robust to the load conditions and number of broken bars; even for challenging scenarios.

  Ês is the (m × d) matrix containing the d eigenvectors associated with the d largest eigenvalues, rearranged in the diagonal matrix Λs whereas Ên is the (m × md) matrix containing the remaining (md) eigenvectors associated with the (md) remaining eigenvalues, rearranged in the diagonal matrix Λn . Obviously, the columns of Ês span a signal subspace estimate.

  and b = [ b0 . . . bdu ] T

[

  S u W T u -I T ] T belongs to the kernel (or identically the null space) of [ Âu Y r ]. We known, by definition, the rank of [ Âu Y r ] is d u column-wise whereas the dimension of this matrix is (m × Nm + 1). Finding [S u W T u -I T ] T is equivalent therefore to seeking for the best d u -rank approximation of [ Âu Y r ]. Owing to the Eckart -Young-Mirsky theorem, the best d u -rank approximation is obtained by a truncated Singular Value Decomposition (SVD) of the matrix [ Âu Y r ]. We have then the decomposition

Figure 1 :

 1 Figure 1: From left to right, healthy rotor and 2 Broken Bars, snapshot of the experimental test Bench

  ) we have first and foremost Hilbert transform the stator current signal that was measured. We first show and remind the difficulties of periodogrambased analysis. To illustrate the concept, we use the signal of a healthy machine under nominal level of load and only 10s of acquisition with sampling frequency equal to 10 KHz. We invite the reader to refer to Fig.2for visual results. As we could expect, the 50 Hz driven frequency is by far the most powerful component inside the spectrum with about 30 dB. From Fig.2-(a) we can clearly see some of the harmonics of the supply frequency, say 100 Hz and 150 Hz with odd harmonics prominent; we also see approximate frequencies of 25,67 Hz and 74,43 Hz. These 2 frequencies are extremely relevant since they correspond to the 50 Hz amplitude modulated by the rotation motion. The slip can be estimated by the side components 52,53 and 47,55 Hz given by Fig.3-(b), which illustrates the proposed method when we have lows-pass filtered the signal and downsampled it by a factor of 10. Owing to the relation (1 + 2g)f r = f sb we can find approximately g = 0, 025 and next deduce the rotation speed of the rotor by 25-25*0.025=24,37 Hz. Then, the 2 frequencies figured out here are[START_REF] Bellini | Advances in Diagnostic Techniques for Induction Machines[END_REF]62 and 74,37 and 50+24,[START_REF] Bouleux | An Optimal Prior knowledge-based DOA Estimation Method[END_REF]. Let us pause and note that the interpretation of these frequencies is only made possible by the use of Fig.3-(b), which clearly gives evidence to the sidebands components. Neither Fig.2-(a) nor Fig.2-(b) have sufficiently discriminated the frequencies to have the information of the sidebands; even if Fig.2-(b) uses the Welch optimal estimator.

Figure 2 :Figure 3 :

 23 Figure 2: 10s of acquisition for a healthy machine, (a) PSD, (b) Welch Estimator with optimal overlapping windows

Figure 4 :

 4 Figure 4: PSD of current signal Low pass Filtered for a healthy machine, (a) without downsampling, (b) with downsampling equivalent to a ZOOMFFT method.

Figure 5 :

 5 Figure 5: 800ms of acquisition, current signal Low pass Filtered and downsampled by a factor of 10 for nominal level of load, 1 broken Bar, (a) PSD (b) Estimated Power Spectrum by the proposed method.

Figure 6 :

 6 Figure 6: 800ms of acquisition, current signal Low pass Filtered and downsampled by a factor of 10 for low level of load, 1 broken Bar, (a) PSD (b) Estimated Power Spectrum by the proposed method.

Figure 7 :

 7 Figure 7: 1s of acquisition, current signal Low pass Filtered and downsampled by a factor of 10 for nominal level of load, 2 broken Bars, (a) PSD (b) Estimated Power Spectrum by the proposed method.

  for 2 Broken Bars -case of 1000 samples (1s of acquisition)

Figure 8 :

 8 Figure 8: 1s of acquisition, current signal Low pass Filtered and downsampled by a factor of 10 for low level of load, 2 broken Bars, (a) PSD (b) Estimated Power Spectrum by the proposed method.

Table 1 :

 1

	State	f r (Hz)	k	2gf r (Hz) I (×10 -4 )
	Healthy	49.98	±1	0.6	3.12
	1 Broken Bar	50.31	±1	1.2	796
	2 Broken Bars	49.82	±{1, 2, 3}	1.4	4147

Table 2 :

 2

  8), we have reported the indicator values inside respectively table 2 and table 3, for nominal level of load and low level of load. We did not dispose however of the healthy signal acquisition for the low level of load scenario which explains why the indicator is absent in table 3.

	State	f r (Hz)	k	|2kg|f r (Hz) (left) (right)	I (×10 -6 )
	1 Broken Bar	50	±1	2.68	2.7	245
	2 Broken Bars	50.01	±{1, 2}	2.48 7.44	2.5 7.5	1500

Table 3 :

 3


Nikolaos Kekatos
email: nikolaos.kekatos@univ-grenoble-alpes.fr

Marcelo Forets
email: marcelo.forets-irurtia@univ-grenoble-alpes.fr

Goran Frehse
email: goran.frehse@univ-grenoble-alpes.fr

Modeling the Wind Turbine Benchmark with PWA Hybrid Automata (Experience Report)

The wind turbine benchmark is part of the ARCH benchmark repository. It entails closedloop requirements and encompasses nonlinear and hybrid dynamics. Owing its origin to industry based applications, the benchmark modeling is done with MATLAB/Simulink. Formal verification tools, however, do not operate on simulation models but on formal models, such as hybrid automata. Particularly efficient verification algorithms are known for systems with Piecewise Affine (PWA) dynamics. In this vein, we construct a PWA model of the wind turbine in the SX format, which formally describes a network of hybrid automata and can be used by several reachability tools. The model transformation follows a four-step approach with the aim of (i) adapting the Simulink model to obtain a verification model (ii) translating the Simulink blocks to equivalent blocks in SX format, (iii) conducting compositional, syntactic hybridization to obtain a PWA approximation of the dynamics of the nonlinear blocks, and (iv) performing model validation. We also report some preliminary experiments on the subsystems (network components) of the wind turbine that we conducted with SpaceEx.

Introduction

The benchmark was introduced by Simone Schuler, Fabiano Daher Adegas, and Adolfo Anta from General Electrics [START_REF] Schuler | Hybrid modelling of a wind turbine (benchmark proposal)[END_REF] in the ARCH workshop. The authors presented a simplified nonlinear model of a wind turbine equipped with switching controllers. The composition of the wind turbine and the controllers results in a hybrid system with nonlinear dynamics. Given that wind turbine systems form one of the fastest-growing industries in renewable energy worldwide, the authors consider them as a promising area for the application of verification tools. The accompanied control objectives including safety and performance constraints establish a pertinent benchmark for verification of requirements in hybrid systems [START_REF] Schuler | Hybrid modelling of a wind turbine (benchmark proposal)[END_REF].

The main issue that arises is that formal verification tools do not operate on simulation models but on formal models. As such, the MATLAB/Simulink model has to be transformed in order for the requirements of the wind turbine benchmark to be verified. This paper describes the process of constructing a Piecewise Affine (PWA) model expressed in SX format (SpaceEx modeling language) [START_REF] Cotton | The SpaceEx modeling language[END_REF] from the Simulink model and presents preliminary verification results undertaken with the SpaceEx tool [START_REF] Frehse | Spaceex: Scalable verification of hybrid systems[END_REF]. The model transformation includes four main steps. We first change the Simulink model to comply with the verification standards (adding nondeterminism, deleting unnecessary blocks). Then, we employ the SL2SX tool, [START_REF] Minopoli | SL2SX translator: From Simulink to SpaceEx models[END_REF] to assist with the translation and the process of building the SpaceEx model. Afterwards, we conduct compositional syntactic hybridization [START_REF] Kekatos | Constructing verification models of nonlinear Simulink systems via syntactic hybridization[END_REF] to obtain a PWA approximation of the dynamics of the nonlinear blocks and finally we perform model validation to check that the individual approximations (base components) are correct and non-blocking.

Much work has been done towards the verification of Simulink models [START_REF] Alur | Symbolic analysis for improving simulation coverage of Simulink/Stateflow models[END_REF] and several tools have been designed to facilitate this task [START_REF] Zhan | Formal Verification of Simulink/Stateflow Diagrams: A Deductive Approach[END_REF]. There are two main groups of approaches. The first one refers to verification by simulation methodologies that can be applied directly on the Simulink model and do not require any model transformation. Donzé [START_REF] Donzé | Breach, a toolbox for verification and parameter synthesis of hybrid systems[END_REF] presented a MAT-LAB/Simulink based tool, Breach, which performs simulation-based verification (approximate reachability analysis) and conducts efficient signal monitoring of properties and requirements. Another relevant MATLAB/Simulink tool is S-Taliro [START_REF] Annpureddy | S-taliro: A tool for temporal logic falsification for hybrid systems[END_REF]. It conducts fast and efficient simulations, with a formal focus on falsification. These two approaches have the drawback that the set of initial states must be sampled and are restricted to models with low-dimensional initial states. A tool that circumvents this pitfall is HySon [START_REF] Bouissou | HySon: Set-based simulation of hybrid systems[END_REF]. It performs set-based simulation on Simulink models, and approximates the set of all possible executions. However, it is not publicly available and, to the best of our knowledge, it has been only tested in cases of relatively small size. The second group of approaches relies on model transformation and requires the construction of a verification model from its Simulink counterpart. Once the verification model is designed, different formal verification tools could be applied. Filipovikj et al. [START_REF] Filipovikj | Simulink to UPPAAL statistical model checker: Analyzing automotive industrial systems[END_REF] transformed Simulink models into the input language of UPPAAL and used the UPPAAL Statistical Model checker to verify properties of automotive systems. The Simulink-based tool C2E2 [START_REF] Fan | Automatic reachability analysis for nonlinear hybrid models with C2E2[END_REF] transforms Stateflow models to polynomial hybrid automata, saves them in the HyXML format, and then conducts simulation-based verification.

The rest of this paper is organized as follows. Section 2 reviews the existing tools that are utilized. In Section 3, we briefly present the benchmark model. In Section 4, we describe the model transformation process. We report some preliminary results in Section 5. We draw conclusions and describe our future work in Section 6. The model and configuration files for running the case study in SpaceEx are available as an attachment.

Background

In this section, we briefly present the three tools that are used for the model transformation, namely, Simulink, SpaceEx and SL2SX, and describe the hybridization technique that we apply.

Simulink

Simulink [START_REF]MATLAB 9.0 and Simulink 8[END_REF] is a graphical programming environment for modeling, simulating and analyzing dynamical systems. It includes a set of block libraries and is commonly used in automatic control and signal processing for multidomain simulation and Model-Based Design. A Simulink model consists of a set of input and outputs variables, blocks and connections.

SpaceEx

SpaceEx [START_REF] Frehse | Spaceex: Scalable verification of hybrid systems[END_REF] is a verification platform for hybrid systems. Its main functionality is to verify that a given mathematical model of a hybrid system satisfies desired safety properties. SpaceEx offers three different algorithms for set-based reachability analysis, namely the PHAVer, the LGG and the STC, as well as a simulation option. The PHAVer (Polyhedral Hybrid Automaton Verifyer) algorithm is applicable to linear hybrid automata, i.e. hybrid systems with piecewise constant bounds on the derivatives, and produces exact reachibility results. The LGG Support Function algorithm implements a variant of the Le Guernic-Girard (LGG) algorithm [START_REF] Le | Reachability analysis of linear systems using support functions[END_REF]. The STC (Space-Time approximation with Clustering) algorithm relies on LGG but it computes fewer convex sets for a given accuracy and yields more precise results for discrete transitions. These algorithms over-approximate the reachable sets and can be applied to hybrid systems with piecewise affine dynamics and non-deterministic inputs.

The models for SpaceEx are specified in the SX modeling language [START_REF] Cotton | The SpaceEx modeling language[END_REF]. The format is similar to the standard hybrid automata, syntactically extended with hierarchy and templates. An SX model consists of base and network components as well as binds. Base components correspond to single hybrid automata, whereas network components correspond to the parallel composition of several hybrid automata. A bind instantiates a component inside a network (similar to a template instantiation), possibly remapping its variables to variables of the network. The parameters of a component can be variables (changing over time), constants or synchronization labels (every transition is associated with a label). The variables can be real-valued or boolean (depending on their type), local or global (depending on their usage outside the base component) and controllable or uncontrollable (typically, distinguishing between variables that appear in differential and algebraic equations). More information can be found at [START_REF] Donzé | Modular, hierarchical models of control systems in SpaceEx[END_REF].

SpaceEx requires two input files to run, a model file of the hybrid automata in SX format and a configuration file that specifies the analysis options and preferences. In the latter file, the user can define, among else, the initial states, output formats, number of iterations, e.g. searching for a fixed point, time horizon, forbidden states. An informal introduction to SpaceEx and its modeling paradigm is available at [START_REF] Frehse | Introduction to SpaceEx v0.8[END_REF].

Simulink to SpaceEx Translator (SL2SX)

SL2SX [START_REF] Minopoli | SL2SX translator: From Simulink to SpaceEx models[END_REF] is a semi-automated tool that undertakes the translation of Simulink models into SpaceEx models. The translator accepts a Simulink model that is saved in XML format and generates as an output a network of hybrid automata in SX format. The translation preserves most of the structural aspects of the Simulink diagram, such as the names, hierarchy and graphical positions. Blocks that are not supported are translated into "empty" components, which must then be completed by the user.

Compositional Syntactic Hybridization

Hybridization is an abstraction method to obtain PWA over-approximations of nonlinear dynamics [START_REF] Asarin | Hybridization methods for the analysis of nonlinear systems[END_REF]. It consists of partitioning the state-space into a set of domains, and for each domain, approximating the nonlinear dynamics by simpler ones plus nondeterministic inputs to account for the abstraction error. Existing hybridization methods operate on the composed (flattened) system. Compositional syntactic hybridization, a recently proposed method [START_REF] Kekatos | Constructing verification models of nonlinear Simulink systems via syntactic hybridization[END_REF], performs the PWA over-approximation in a compositional manner and takes advantage of the on-the-fly composition of hybrid systems that is supported by the SpaceEx platform. Three main steps are involved: syntactic decomposition, replacing the nonlinear ODE by a linear ODE with nonlinear algebraic equations; hybridization, constructing a PWA approximation for each algebraic equation and providing a sound over-approximation of the original system by adding an error term; and finally composition of the hybrid automata, and elimination of the algebraic equations.

Benchmark Model

The wind turbine modeling is done with MATLAB and Simulink [START_REF] Schuler | Hybrid modelling of a wind turbine (benchmark proposal)[END_REF]. Figure 1 shows the top-level model of the Simulink block diagram and Figure 2 depicts the physical part (plant) of the wind turbine. The wind turbine dynamics are highly nonlinear functions of the operating point defined by the rotor speed, wind speed and blade pitch angle.

The plant consists of three subsystems: a servo-elastic, an aeroelastic and a pitch-actuator. The servo-elastic subsystem describes the tower fore-aft dynamics and the rotor dynamics. It is interconnected with the aeroelastic subsystem through signals that correspond to the aerodynamic torque and thrust. While the servo-elastic subsystem consists of linear operators, the aeroelastic subsystem has several nonlinear blocks. There are two Square functions, two Products, two Divisions and two fourth-order Polynomials. The polynomials correspond to the aerodynamic torque coefficient (cP ) and thrust coefficient (cT ), respectively. They are computed through a regression model. The pitch actuator dynamics can be expressed by a second-order lag, a first-order lag or a time-delay. The pitch actuator subsystem has two input signals, the demanded pitch angle and the pitch actuator type and outputs the actual pitch angle to the aeroelastic subsystem. Apart from the linear operators, it contains a Multiport Switch, three Compare to constant blocks, three Enable blocks and three Enabled Subsystems.

As for the control part, there is a generator-torque controller and a collective blade-pitch controller. The blade pitch controller is a gain-scheduled PID whose objective is to minimize the speed error between the filtered generator speed and the rated generator speed. It incorporates an anti-windup scheme to prevent integration wind-up when the actuator is saturated. The nonlinear blocks correspond to a Division and a Product. The generator-torque controller aims to maximize the extracted maximum power from the wind by tracking the optimal tip-speed ratio λ opt . It is a hybrid controller with 5 locations and 2 input signals, namely the filtered generator speed and the pitch angle. It is written as an Embedded MATLAB Function and it includes two nonlinear operators, a Square and a Division.

Note that the benchmark model is designed from an industrial viewpoint and its verification task is considered to be of high difficulty. That ensues from the existence of several nonlinearities (over 10 distinct ones) and blocks that cannot be expressed by linear or hybrid dynamics. As far as the model dimension is concerned, the controllers introduce two state variables and the plant introduces 3 -5 states (depending on the pitch actuation). As such, the total number of state variables of the Simulink (closed-loop) model varies from 5 to 7. This number is further increased, if we consider the rate limiters and memory blocks. The SpaceEx model adds an extra variable to capture the time evolution. Another point is that the Simulink model is more detailed than the mathematical model described in the benchmark paper [START_REF] Schuler | Hybrid modelling of a wind turbine (benchmark proposal)[END_REF] (e.g. rate limiters). There are two minor points regarding the benchmark model that are worth mentioning; some MATLAB paths need to be modified in order to be platform/user independent and there is a call to a toolbox (WAFO) that is not included in the benchmark files.

Model Transformation

This section presents the model transformation process, covering the model adaptation, translation, hybridization and validation steps.

Model Adaptation

To apply formal verification techniques, it is necessary to obtain a verification model from the original simulation model. Such a transformation should consider the inherent differences between these two model classes. Typically, a simulation model includes details that should be obscured from the verification model. A verification model could also be enriched with nondeterminism so as to check the behavior of the system under uncertain or varying parameters, disturbances or user inputs.

An illustrative example is the Simulink block corresponding to the wind speed profile. The wind profile is an input to the model that is read from a data structure in the MATLAB workspace. In principle, such a profile could be translated into a hybrid automaton, but this would not be very efficient in the analysis. As a matter of fact, there are 33 possible profiles of user-defined wind speed signals in the MATLAB aeromaps.mat file. However, even if all these cases are tested, there is still no guarantee that the closed-loop system is going to operate correctly. Actually, there might be a different wind profile that could yield undesirable behaviors. On the contrary, under the SX format, this matter could be simply resolved by adding a non-deterministic input signal that covers the minimum and maximum wind speeds. It is possible to impose bounds on the rate of change of the wind speed.

It is also necessary to preprocess and further modify the Simulink model. In fact, all the Scopes, Mux, Demux, Enabled Subsystems, Manual Switches, Save to workspace have to be deleted/replaced as their corresponding actions are defined in a different way with SX format and SpaceEx1 . Another modification concerns the three available pitch actuators that exist in the Simulink model. As the pitch actuator dynamics cannot be neglected for large wind turbines, we consider the most critical and representative case, i.e. the second-order lag.

Translation

The second step of the model transformation is to translate the Simulink model into an equivalent SpaceEx model; equivalent in the sense that the model remains the same while the syntax changes. For this task, we employ the SL2SX translator. Having the mechanical aspects of the model translation carried out by a tool significantly reduces errors. The translation takes a few seconds and produces an SX file with 18 Network and 89 Base Components. However, not all Simulink blocks can be exactly translated, as they cannot be expressed as hybrid automata. The Simulink blocks that are automatically translated into base components are the following: Add, Subtract, Divide, Multiply, Constant, Gain, Saturation, Integrator, Subsystem, Inport, Outport.

Finally, it is possible to make the SX model more compact by reusing base components and deleting the ones that are duplicate or serve the same role. With this configuration, the total number of base components can be reduced approximately up to 30 distinct ones2 . Figure 3 and Figure 4, respectively, present the top level block (network component) and the plant subsystem of the wind turbine model in SX format, as shown in the Model Editor of SpaceEx.

Implementation Issues After modifying the Simulink blocks, we have to perform further preprocessing in order to be able to use the SL2SX translator. Special attention should be given to the names and the connections. As of version 1.0.1 of SL2SX, the block names should not extend to multiple lines, as the respective blocks will not be parsed by the translator. There may occur errors with models that are saved in old Simulink versions. Also, after constructing the SX Model with the SL2SX tool, there is a need for post-processing. In particular, the translator produces the expressions "; &gt" and "; &lt" instead of the correct inequality symbols "<" and ">", respectively3 .

Hybridization

The third step is to generate PWA approximations and describe them in SX format for all these Simulink blocks that cannot be handled automatically, either because no exact translation is available (e.g. nonlinearities) or because translation cannot be applied (e.g. Embedded MATLAB Function).

To obtain PWA approximations of nonlinear dynamics, we use an abstraction method we call syntactic hybridization. This method essentially partitions the state-space into a set of domains, and for each domain, it approximates the nonlinear dynamics by simpler ones with added nondeterministic inputs to account for the abstraction error. However, instead of operating on the composed (flattened) system, we decompose the original dynamics and carry out the state-space partitioning and PWA approximation on the components. In this way, we can avoid having intractably large models and the explosion in the number of partitions can be mitigated. In essence, we break down the nonlinear blocks into components that have a small number of input variables and interconnect them together.

A practical intermediate step is to get bounds on the behavior (min, max) of the input signals of the nonlinear Simulink blocks. The shorter the ranges of the signals, the smaller the number of locations that is required by the PWA abstraction, given a desired error bound. There are different ways to estimate them. In this paper, we run Simulink simulations for different scenarios and initial conditions to obtain a (not necessarily conservative) estimation of the signal range. Note that the signal range serves only as an indication for the hybridization step. The approximation is equipped with out-of-range (error) states. So, if the range is found to be insufficient during reachability computations, it is revised (enlarged).

For the wind turbine benchmark, the constructed base components have only one or two input signals and one output signal, the state space is partitioned into a set of domains, described by hyper-rectangles, the approximations are linearized (first-order Taylor) around the center of each domain, and the abstraction (linearization) error is computed by the evaluation of the Lagrange remainder. The maximum error is computed for each location (box). We use a combination of interval arithmetic and global optimization to bound it [START_REF] Althoff | Reachability analysis of nonlinear systems with uncertain parameters using conservative linearization[END_REF]. All the computations (linearization domains, quantization parameters, operating points, PWA approximations, errors) are conducted with MATLAB. We also have a script to facilitate the construction of PWA approximations and transform them into equivalent base components models, described in SX format. Then, we integrate all these components into a single SX file.

Note that the approximations can be conservative or non-conservative. The error is computed for each location of the base components and is added as a non-deterministic input in the corresponding invariants. However, it is possible to disregard the errors either for debugging purposes or for the purpose of obtaining a deterministic model. Table 1 presents all the Simulink blocks that need to be approximated along with the approximation method, the maximum induced errors and the number of locations of each approximation. Let us underscore that the error corresponds to the maximum error that appears in one of the locations of the PWA approximation. The individual errors in the other locations can be significantly smaller.

The resulting model has only 72 locations in all components combined. To carry out reachability analysis, SpaceEx performs on-the-fly composition of the hybrid automata and instantiates only the reachable parts of the state-space. The upper bound of the number of locations of the composed model is 16 millions.

Model Validation

The fourth step of the model transformation concerns the model validation. We propose an empirical method to evaluate whether the generated base components are correct, insofar as they yield satisfactory behaviors and do not introduce any deadlocks. In this step, we are essentially testing the implementation. In this vein, we have created an SX signal library with predefined components that describe trigonometric (sine, cosine), step and ramp functions. We have added a monitor to visualize the output 4 . Once we create a new approximation, we can integrate it along with the SX library and obtain a tester module (consisting of base and network components).

Our objective is twofold: on the one hand, to check that the approximation is non-blocking and on the other hand to check that it is indeed an over-approximation. For the first objective, the input signals are selected in a way that compels the PWA automaton (introduced by the hybridization process) to necessarily visit all of its locations. As soon as the input signals exceed the allowed operating range, the automaton goes out-of-bounds and the analysis terminates. That means that a fixed point has been reached. For the second objective, we perform reachability analysis on the tester module (containing the approximation and the input signals) and then compare the results with random simulations. In this way, we check whether the simulations are included in the reachable sets. As an illustration, let us consider a Simulink block of the aero-elastic subsystem of the wind turbine (plant). Figure 5 shows the entire aero-elastic subsystem. The block that we focus on is the Divide block and it is highlighted in red. This block corresponds to the division operator and computes the division between two input signals over time. From a physical point of view, it computes the tip-speed ratio, which is a dimensionless variable denoted by lambda (λ). This block, being a nonlinear one, is replaced by a PWA over-approximation with 4 locations. Figure 6 shows the approximation as a SpaceEx base component. The variables In1 and In2 relate to the inputs, Out1 is the PWA approximation of the division operator In1/In2 and corresponds to λ. The approximation error w1 is added in the form of a non-deterministic disturbance. Note that all the variables should be uncontrollable, and the error term should be included in the invariant, declared as a local variable; see Section 2.2 for more information about SpaceEx modeling language.

After incorporating the base component of the Divide block with our input library, we want to check the previously mentioned objectives (no deadlocks and over-approximation). For this reason, we consider that the inputs are ramp signals with varying initial conditions. The first input is described by dIn1 dt = 1 with initial conditions in 69 < In1 < 71 and the second input by dIn2 dt = 1 with initial conditions in 6.9 < In2 < 7.1. It is necessary to add a monitor in order to visualize the output, as we intend to use the STC algorithm. Then, we conduct reachability analysis, compare the results with random simulations of the original nonlinear function and observe that for the considered scenario SpaceEx yields an over-approximation. The tester module is shown in Figure 7a, the reachable sets and the simulation runs are shown in 7b. The same approach (in a hierarchical way) is utilized in order to check that compositions of multiple base components, network components or larger subsystems operate correctly.

Preliminary Results

In this section, we present validation runs for some of the subsystems (network components). We start with the pitch-actuator subsystem of the wind turbine plant. The input of this block is the commanded pitch angle and the output is the actual pitch angle. Considering the STC scenario, a global time horizon of 10s, and a constant input signal, we get a fixed point. Figure 8 depicts the reachable sets for different initial conditions. The input signals are shown in red and the output signals in blue.

Next, we consider the servo-elastic subsystem of the wind turbine plant. Given constant inputs, a global time horizon of 20s and the STC scenario, SpaceEx finds a fixed point. Figure 9 depicts the reachable sets for initial conditions: 0 ≤ Omega ≤ 0.1 and 0 ≤ xT dot ≤ 0.1. Finally, we analyze the torque controller. The approximation consists of 28 locations, has two input signals (omega d, theta) and one output (torque). We consider the input signals of Figure 11 (a ramp and a sine function), as they cover, when combined, the entire 2D operating range (from minimum to maximum allowed values). Indeed, with this configuration, SpaceEx visits all the controller locations and finds a fixed point after 16 seconds. The reachable sets are computed with the PHAVer scenario. The reachable sets of the torque controller are displayed in Figure 10.

Conclusions

Model transformation plays an important role in bridging the gap between industrially relevant models and verification tools. This work aims to assist the application of hybrid system reachability tools to industrial-sized models described by MATLAB/Simulink. The existence of a translator from Simulink to SpaceEx already facilitates the use of Simulink models within SpaceEx environment. However, the large variety of Simulink blocks and their diverse features pose difficulties in mapping Simulink subsystems to SpaceEx components. There are several Simulink blocks that cannot be described by hybrid automata with PWA dynamics and therefore cannot be processed by SpaceEx. In practice, most of these blocks correspond to nonlinear functions and operations.

The wind turbine benchmark constitutes a relevant example, insofar as it represents an industrial case study, while including many nonlinear blocks. Through syntactic hybridization, it is possible to over-approximate the dynamical behavior of these blocks and seamlessly integrate them within SpaceEx components. Note that the resulting model of hybrid automata is expressed in the general SX format. As such, it can be fed directly into SpaceEx platform, or translated into formats that comply with other verification tools using the HyST [START_REF] Bak | HYST: a source transformation and translation tool for hybrid automaton models[END_REF] tool.

Our future work is targeted towards verifying the requirements of the wind turbine benchmark. As the proposed approximation contains around 16 million locations, there is a shortcoming concerning the computational costs and execution time. Currently, SpaceEx identifies the initial locations of the approximation through enumeration, which does not scale. A potential way to tackle this issue is to use compositional reasoning to identify the initial locations and instantiate as few locations as possible during the analysis.

Figure 1 :

 1 Figure 1: Wind turbine -Simulink model -Top Level [18].

Figure 2 :

 2 Figure 2: Wind turbine -Simulink model -Plant [18].

Figure 3 :

 3 Figure 3: Wind turbine -SpaceEx model -Top Level.

Figure 4 :

 4 Figure 4: Wind turbine -SpaceEx model -Plant.

Figure 5 :

 5 Figure 5: Aero-elastic subsystem -Simulink -Nonlinear Divide block shown in red.

Figure 6 :

 6 Figure 6: Divide block -PWA approximation in SpaceEx.

  (a) Tester -Division λ. (b) Reachable sets and Simulations.

Figure 7 :

 7 Figure 7: SpaceEx Base Component -Division λ -Model Validation. The reachable sets of the approximate SpaceEx block are indicated with blue, the simulation runs of the nonlinear function are denoted in red. The simulation runs are contained in the reachable sets.

  (a) Initial Conditions -0 ≤ θ ≤ 0.15. (b) Initial Conditions -0.18 ≤ θ ≤ 0.22.

Figure 8 :

 8 Figure 8: Pitch actuator -SpaceEx Network Component -Reachable Sets (blue: Output, red: Input)

Figure 9 :

 9 Figure 9: Servo-elastic -SpaceEx Network Component -Reachable Sets

Figure 10 :

 10 Figure 10: Torque Controller -SpaceEx Network Component -Output.

(a) Input 1 -

 1 Omega d. (b) Input 2 -Theta.

Figure 11 :

 11 Figure 11: Torque Controller -SpaceEx Network Component -Inputs.

Table 1 :

 1 SpaceEx base components -Approximationsa for visualization reasons an extra location should be added.

	No Block (Name)	Approx.	# Loc. Error Bounds Remarks
	1	Product (Anti windup) syntactic	4	2	x • y
	2	Division (C p /λ)	syntactic	10	0.027613	(x/y)
	3	Division (GS factor )	syntactic	2	0.020111	1 1+x
	4	Division (λ)	syntactic	4	0.6983	x/y
	5 6	Product (C T • v 2 ref ) Product (C p /λ • v 2 ref )	syntactic syntactic	4 4	25.5781 3.23437	x 2 • y x 2 • y
	7	Polynomial (aeromaps) syntactic	6	3.39472	4th-order
	8	Polynomial (aeromaps) syntactic	6	12.2809	4th-order
	9	Embedded MATLAB	syntactic	28	10.48584	x 2 and 1/x
	10 Saturation	substituted	3	0	exact a
	11 Read from workspace	substituted	1	-	range b
	12 Mux, Scope	unnecessary -	-	-
	13 Save to workspace	unnecessary -	-	-
	14 Enabled Subsystem	unnecessary -	-	-
	15 Multiport Switch	unnecessary -	-	-
	16 Compare to Constant	unnecessary -	-	-
	17 Manual Switch	unnecessary -	-	-
	18 Rate Limiter	ignored c	-	-	DAE
	19 Memory	ignored c	-	-	continuous delay

b replaced by a nondeterministic input. c to mitigate state explosion, this block has not been considered.

Currently, the preprocessing is done manually. However, in principle, this step could be automated. Deleting unnecessary blocks and checking naming issues should be a straightforward task.

This process is currently manual but it could be automated; either as a functionality of SL2SX or as a post-processing step.

This a known issue with the way that SAX parser (which is called by SL2SX) handles special characters.

Note that STC and LGG algorithms in SpaceEx do not output algebraic variables, so visualizing is possible by mapping the algebraic variable to the solution of an ordinary differential equation.

Acknowledgments

The authors would like to thank Simone Schuler for valuable discussions on wind turbine models.

The authors gratefully acknowledge financial support by the European Commission project UnCoVerCPS under grant number 643921 and by the Metro Grenoble through the project NANO2017.