Pierre-Louis Chiambaretto

Miguel Charlotte

Joseph Morlier

Philippe Villedieu

Yves Gourinat

Surrogate Granular Materials for Modal Test of Fluid Filled Tanks

Keywords:

des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

INTRODUCTION

Cryogenic stages are commonly used for launchers due to their efficiency. However, several technical issues make their development tough. One of these challenges is to determine the dynamic responses of tanks entirely or partially filled with cryogenic fluids. Indeed, it is dangerous and expensive to perform vibration tests on tanks filled with these fluids. Therefore, for some fluids such as liquid oxygen, surrogate fluids with close thermodynamic constants can be used. As explained by NASA [START_REF] Tuma | Ensuring Safe Exploration: Ares Launch Vehicle Integrated Vehicle Ground Vibration Testing[END_REF] , this approach is also used for liquid hydrogen, by using water even if the density ratio of these fluids is around 14. Looking for alternatives, this innovative work deals with the possibility to use pre-stressed granular materials to fill the tank in order to obtain experimentally the same modal response as a tank filled with liquid hydrogen for low frequencies.

A methodology based on a semi-analytical resolution to choose the granular material and the pre-development of an experimental test bench is presented: After the justification of 2D hypothesis, the aim will be to validate the modal equivalence between fluid filled tanks and grains filled tanks. In order to achieve this, another model based on a homogenization of the granular material is used as sketched in figure 1 . After explaining the hypotheses associated with this homogeneous model, eigen-vectors and eigen-frequencies of the proposed mechanical model are determined thanks to the Rayleigh-Ritz method. To obtain the same modal parameters as fluid filled tank, grain material can be chosen. A test bench is under development to confirm this approach and identify some parameters such as the coordination number of the granular material.

beam + fluid beam + membrane beam + grains ? ?

Figure 1: modal equivalence between fluid filled beam and grains filled beam using an homogeneous model for granular material

MODELING

This work focuses on cylindrical tanks such as the Ariane 5 cryogenic main stage tank. The first step is to define a model of tank filled with liquid hydrogen, in order to obtain an expression of the first eigen modes intended for the equivalence with the tank filled with grains.

Fluid-structure interaction

The tank is mainly composed of a long hollow cylinder closed by two spherical caps. Since the thickness is considerably smaller than the other dimensions, consequently the tank can be modeled by a shell. Moreover, the two spherical caps are non-developable surfaces, therefore the rigidity of theses caps is much larger than the cylinder's one. We assume the caps do not participate in the first modes. Then, the model can be reduced to a cylindrical shell clamped. The shell geometry is defined by the radius R and the thickness h and its homogeneous, isotropic material by the Young modulus E1 and the density ρ1. The first modes of this structure are flexural radial-orthoradial modes. Regardless of the shell theory chosen, with the parametrization of the figure 2, using the method of separation of variables and the invariance to rotation, the displacement can be written in the form :

           u(θ, x, t) = i,j uij(θ, x, t) = i,j F u ij (x) cos (iθ) cos (ωt) v(θ, x, t) = i,j vij(θ, x, t) = i,j F v ij (x) sin (iθ) cos (ωt) w(θ, x, t) = i,j wij(θ, x, t) = i,j F w ij (x) cos (iθ) cos (ωt)
.

(1)

The relation between i, j and ω depends on the shell theory used. In order to solve analytically the problem, an assumption of plane motion is used; accordingly u(θ, x, t) v(θ, x, t) and w(θ, x, t). Then, the tank can be studied slice by slice with an amplitude function of x. A tank's slice is modeled by a circular beam. As the radius is much larger than the thickness, the beam theory used is the Euler-Bernoulli's one [START_REF] Love | A treatise on the mathematical theory of elasticity[END_REF] ; Love determined the flexural modes of a circular beam. These modes were obtained by adding to the Euler-Bernoulli equations, a closure relation imposing: θθ = 0 ⇔ ∂v ∂θ -w = 0. This relation imposes that the work done by normal forces (membrane forces) is null. The i th eigen-modes functions and eigen-frequencies equations are:

       vi(θ, t) = Vi(θ) cos (ωt) = Ax sin (iθ) cos (ωt) wi(θ, t) = Wi(θ) cos (ωt) = iAx cos (iθ) cos (ωt) ω = i(i 2 -1) R 2 (i 2 +1) 1/2 E 1 h 12ρ 1 1/2 ∀i ∈ N . (2)
Thanks to these expressions, the generalized mass and stiffness for each mode by unit of length can be obtained as follows:

       K B g |i = E 1 h 3 12R 3 2π 0 d 2 W i dθ 2 + Wi 2 dθ = Ax 2 E 1 h 3 π 12R 3 1 -i 2 2 M B g |i = ρ1hR 2π 0 V 2 i + W 2 i dθ = Ax 2 ρ1hRπ 1 + 1 i 2 . ( 3
)
In this work, the tank is assumed totally filled, so there is no sloshing effect. To determine the importance of acoustic effects, orders of magnitude are used. The celerity of sound c is about 1000 m/s, the first frequency f of our system is about 10 Hz and the characteristic 2R is about 1 m. The Helmholtz number is defined for example in [START_REF] Sigrist | Interaction fluide-structure: Analyse vibratoire par él éments finis[END_REF] as :

He = 2πf R c 2 ≈ 10 -2 .
As a consequence, the hypothesis H1 is made :

• H1 : Acoustic effects are negligible.

The Helmholtz equations simplify to Laplace equation in pressure [START_REF] Sigrist | Interaction fluide-structure: Analyse vibratoire par él éments finis[END_REF] :

∆(p) = 0 . (4)
To model the entire system, equations of elastic beams and equation 4 are used. The interaction is defined on the boundary of the liquid domain by :

σ( n) = -p n --→ grad(p). n = ρF ω 2 w (5)
n representing the normal vector from the fluid to the beam and ρF the density of the fluid. As shown for example by Axisa in [START_REF] Axisa | Modelling of Mechanical Systems: Fluid-Structure Interaction: Fluid-Structure Interaction[END_REF] for each mode, the mode shapes remain the same than for the beam without fluid, and the frequency decreases due to the inertia added by fluid. This phenomenon can be decreased by an added mass due to the fluid on each mode. In the case of a ring, this added mass by unit of length is defined by : M F a |i,j = F w ij (x)

2 ρ F πR 2 i
. Here F w ij (x) denotes the amplitude of the beam radial motion for the slice. In this 2D approach, it is a constant. Eigen-frequencies, can be obtained thanks to Rayleigh ratio with

M F a |i = M F a |i,j with F w ij (x) = 1 : fi = 1 2π K B g |i M B g |i + M F a |i . ( 6
)
Thanks to this analysis, it is possible to sum up on modal objectives for the tank filled with grains as follows:

• mode shapes of the tank filled with grains must be the same as the empty tank's one;

• eigen-frequencies must follow relation 6, this paper focuses only on the first mode.

Homogenization of elastic granular materials

A granular material is composed of a lot of grains. The understanding of its different behaviours and the transitions between them, for such materials, constitutes an important research field. In this study the granular material is in its "static material" state: meaning the grains cannot flow. As explained for example by Nedderman [START_REF] Nedderman | Statics and kinematics of granular materials[END_REF] , this state can exist if the material is pre-compressed and the amplitude of pressure oscillations is low enough to prevent the formation of slip planes. According to the Mohr-Coulomb theory [START_REF] Nedderman | Statics and kinematics of granular materials[END_REF] , this failure mode happens when :

1 -sin(δ) 1 + sin(δ) σ1 > σ2 > 1 + sin(δ) 1 -sin(δ) σ1 . (7)
Here σ1 and σ2 are the principal stresses (the third principal stress being null in this 2D approach, the same kind of relation exists in 3D). δ is the effective angle of internal friction computed thanks to the geometry of grains and their assembly and to the friction coefficient between grains.

Given the foregoing hypothesis on principal stresses, the granular material can be described like a non-linear elastic solid. Indeed, even if the grain material is linearly elastic, according to Hertz theory [START_REF] Hertz | Über die Ber ührung fester elastischer K örper[END_REF] the stiffness is non-linear due to the increase of contact area. If the grain's sizes is very small in front of the tank size, an homogeneous theory could be used. In this preliminary work, Boussinesq theory [START_REF] Boussinesq | Essai th éorique sur l' équilibre des massifs pulv érulents: compar é à celui de massifs solides, et sur la pouss ée des terres sans coh ésion[END_REF] is used. The behaviour law is Hooke's law with material parameters depending on the pre-compression of grains : σij = K(∆)∆δij + 2G(∆) ij with ∆ the volumic variation : ∆ = -δV V =kk . Material coefficients can be approximated combining both Hertz contact theory in 2D and Duffy-Mindlin model [START_REF] Duffy | stress-strain relation and vibrations of a granular medium[END_REF] (pressure P is evenly distributed on contacts) as follows : K A1ZE∆ 1/2 A1(ZE) 2/3 P 1/3 and G A2ZE∆ 1/2 A2(ZE) 2/3 P 1/3 . Where Z is the coordination number characterizing the average number of contacts by grain. The Z value depends on the assemblage configuration and the size disparity of grains. For a 2D model, Z ∈ [3, 4] and for a 3D model, Z ∈ [4, 6] as explained by Andreotti and al [START_REF] Andreotti | Les milieux granulaires-Entre fluide et solide: Entre fluide et solide[END_REF] . A1 and A2 are constants depending on grain materials. To sum up, three hypotheses are made :

• H2 : Pre-stressed pressure is high enough to prevent relative motions between grains.

• H3 : Granular material can be described as an elastic material with material coefficient depending on pressure.

• H4 : Pressure variation during vibration tests is low enough to linearize material coefficients.

The equivalent solid is described by its density ρ2, its Lam é parameters λ2 and µ2. Referring to figure 1, the modal equivalence between a "beam+grains" model and a "beam+membrane" model in 2D is now possible. To obtain the modal equivalence between a "beam+membrane" model and the "beam+fluid" model in 2D, it's necessary to determine the mode shapes and eigen frequencies of the "beam+membrane" model. Thus the granular material properties and the pre-stress pressure could be determined to fit with this model's responses to the fluid structure developed in paragraph 2.1

2D Analytical modal determination of a circular membrane surrounded by a beam

The methodology developed to approximate the eigen modes of a circular membrane surrounded by a beam is based on the Rayleigh-Ritz method. A mode shape form is assumed, then by calculating the mass and stiffness operator, eigen-frequencies are determined. To validate the supposed mode shape form, these frequencies are compared to numerical results. With the currently used assumption of separation of spatial variables (r, θ) for mode shapes of disk and using the 2π-periodicity of the structure on θ, the following hypothesis H5 can be made :

• H5 : The radial displacements of the mode shape of a free circular beam, of a circular membrane on its free edge and of a beam filled by this circular membrane are the same.

To validate this assumption, a modal analysis has been made by the finite element method using COMSOL. The figure 3 shows the third mode of a beam, a membrane and a membrane circled by a beam obtained. The deformed shape of the free beam, of the membrane on free edge and of the beam on the coupled system seem to be quit close. Therefore, mode shapes of a free circular beam are the same as for a beam filled by a membrane. This result is coherent with the fluid filled beam. As the mode shapes of the free beam and of the free membrane are barely compatible, for the coupled system these two mode shapes are used for the Rayleigh-Ritz approach with an holonomic condition of same amplitude in radial motion. Mode shapes of the beam have already been computed in section 2.1 as well as the mass and stiffness operators K B g |i M B g |i. In-plane mode shapes of a free edge circular membrane have been less studied than out-plane vibration. In plane mode shapes of a circular membrane with elastic boundary conditions have been analytically obtained by Kim and al. [START_REF] Kim | Exact solutions of in-plane natural vibration of a circular plate with outer edge restrained elastically[END_REF] . By using the separation of variables in polar coordinates, the i th eigen-vector can be written as :

U =   Ur(r)|i cos iθ U θ (r)|i sin iθ 0   , with :          Ur(r)|i = R 2 {A1 Ji+1(α r R ) -Ji-1(α r R ) + A2 Ji+1(cα r R ) + Ji-1(cα r R ) } U θ (r)|i = R 2 {A1 Ji+1(α r R ) + Ji-1(α r R ) + A2 Ji+1(cα r R ) -Ji-1(cα r R ) } c = 2/(1 -ν1) α = ωR ρ(1 -ν 2 2 )/E2 . ( 8
)
where Ji is Bessel function of first kind of order i. The values of α are determined thanks to the dispersion relation :

2Jn+2(α)Jn-2(cα) + 2Jn+2(cα)Jn-2(α) -2(c 2 -1)Jn(α) [Jn+2(cα) + Jn-2(cα)] = 0 . (9)
The amplitudes A1 and A2 are linked by the relation : D22D12A2 = -A1D11D21 with :

         D11 = -λ 2 Ji+2(λ) + Ji-2(λ) -2(c 2 -1)Ji(λ) D22 = -cλ 2 [Ji+2(cλ) + Ji-2(cλ)] D12 = -cλ 2 [Ji+2(cλ) -Ji-2(cλ)] D21 = -λ 2 [Ji+2(λ) -Ji-2(λ)] . (10)
Using classical elasticity theory, the stiffness and mass operators per unit of length can be obtained as follows:

K M g |i = π R 0 1 r [r 2 (U r |i) 2 (λ2 + 2µ2) + 2λ2r(Ur|i + iU θ |i)U r |i + µ2r 2 (U θ |i) 2 -2µ2rU θ |i(iUr|i + U θ |i)+ Ur| 2 i ((2 + i 2 )µ2 + λ2) + 2i(λ2 + 3µ2)Ur|iU θ |i + (i 2 + (λ2 + 2µ2) + µ2))U θ | 2 i ]dr (11)
and

M M g |i = πρ2 R 0 U 2 r + U 2 θ rdr . (12)
To evaluate U r |i and U θ |i recurrence relation on Bessel functions is used :

Ji+1(x) -Ji-1(x) = -2 dJ i (x) dx with x = α or x = cα. U r (r)|i = α 4 {A1 2Ji(α r R ) -Ji-2(α r R ) -Ji+2(α r R ) + cA2 Ji-2(cα r R ) -Ji+2(cα r R ) } U θ (r)|i = α 4 {A1 Ji-2(α r R ) -Ji+2(α r R ) + cA2 2Ji(cα r R ) -Ji-2(cα r R ) -Ji+2(cα r R ) } . ( 13
)
Integration of M M g |i and K M g |i are made numerically. Knowing M B g |i , K B g |i , M M g |i and K M g |i, it is possible using Rayleigh-Ritz theory to determined eigen frequencies of the circular membrane surrounded by the beam:

fi = 1 2π K B g |i + K M g |i M B g |i + M M g |i (14)
For the first mode of ovalization, this analytical frequency is compared to the frequency obtained by the finite element model for different radius and materials of membrane and beam. The maximal error obtained is around 5%. Therefore, this model can be used to find a modal equivalence between the circular membrane surrounded by the beam and the fluid filled beam.

METHODOLOGY OF CHOICE FOR GRANULAR MATERIAL AND FUTURE EXPERIMENTAL TEST BENCH

The objective is to choose a granular material in order to get the same eigen frequencies. Using expression ( 6) and ( 14), this objective can be written as :

K B g |i + K M g |i M B g |i + M M g |i = K B g |i M B g |i + M F a |i ⇔ K M g |i = K B g |i M B g |i + M F a |i M M g |i + K B g |i M B g |i + M F a |i M B g |i -K B g |i . ( 15
)
Fluid acts only as an added mass for the beam whereas the membrane acts as an added mass and an added stiffness. Therefore, the density of the membrane must be greater than the one of the liquid in a proportion depending on the membrane's stiffness as shown in relation (15). To choose the grain material, a density-Young modulus Ashby diagram [START_REF] Ashby | Materials selection in mechanical design[END_REF] , on the figure 4, can be used. This graph represents the main material families, densities and young modulus. For a granular material, the average density ρ2 is linked to the density of the material ρ thanks to the packing fraction c by the relation : ρ2 = cρ. The packing fraction value depends on the organisation of the grains. For a random compact arrangement, the value is closed to 0.7. Young modulus of a granular material E2 is proportional to E 2/3 with E the material Young modulus. Moreover, K M g |i is proportional to E2 and M M g |i is proportional to ρ2. To verify equation (15) with existing materials, Young modulus needs to be quite low. According to Asbhy diagram, foams and elastomer could be interesting. By calculating the mass and stiffness operators of the membrane, for an aluminium beam with a radius of 1 meter, it is possible to determine the pre-stress pressure on grains needed to verify (15). To finalize the choice the pre-stress pressure obtained must be verified at several points during vibration tests:

• Each grain must remain in its elastic domain : σmax < Re ⇒ PMax • Pressure must prevent slipping grain motions : equation 7 ⇒ PMin In order to validate the methodology and the semi-analytical approach, a test bench is under development. Even if the model is mainly 2D, it is quite difficult to realize a 2D experimentation. The model presented on figure 1 of a granular material surrounded by a beam is unstable due to the possibility of out of slice-plane motions. To deal with this issue a solution is to use a long cylindrical tank to be in plane strain state for slices far away from boundaries conditions. The tank is fully filled with spherical grains pressurized by a membrane to get, thanks to a Jassen theory [START_REF] Janssen | Versuche über getreidedruck in silozellen[END_REF] , the pressure needed for the slice of interest. The tank is embedded on a horizontal plate directly connected to the shaker. Accelerometers are settled on the slice of interest. The figure 5 represents a schematic design of the test bench. This work has presented a new approach to perform vibration tests on liquid hydrogen filled tanks. Due to liquid hydrogen's low density, the use of surrogate fluid, as done up to now, presents some important issues. By using a granular material instead of water to fill the tank for vibration tests, it seems possible to reach mode shapes and eigen-frequencies which are unreachable with water.

Prestress accelerometers

After the justification of a 2D approach, the first mode of a cylindrical tank filled with liquid were analytically obtained under the hypothesis H1(incompressible fluid). Once these objective modes are known, an homogeneous approach of granular material based on several hypotheses (mainly H2, H3 and H4) has been developed in order to semi-analytically find the modes of a tank filled with grains using assumption H5. At last, in order to fit the modes obtained for each tank, a methodology for choosing material and pre-stress pressure has been exposed. A test bench, currently under development, has been presented to validate this analytical approach. Theories used in this study are easy-to-apply like the Boussinesq's one. This choice was made to facilitate the analytical approach in order to get orders of magnitude to choose grain properties. The test bench could help us determine the validation range of these theories and to improve modeling if necessary. In future works, we will focus mainly on the impact and the modeling of damping in granular material using experimental results.

Figure 2 :

 2 Figure 2: Cylindrical shell parametrization : ( u, v, w) is the local base of a current point defined by its cylindrical coordinates : (R, θ, x)

Figure 3 :

 3 Figure 3: third mode of a circular beam, a circular membrane and circular membrane surrounded by a beam obtained by COMSOL to validate the hypothesis H5

Figure 4 :

 4 Figure 4: Density-Young modulus Ashby diagram obtained using CES-EduPack.

Figure 5 :

 5 Figure 5: Schematic design of the test bench


Fellow, IEEE Yuri I Abramovich

Senior Member, IEEE Olivier Besson

Ben A Johnson

Conditional Expected Likelihood Technique for Compound Gaussian and Gaussian Distributed Noise Mixtures

Keywords: Direction of arrival estimation, expected likelihood, noise mixture

Expected likelihood (EL) technique for quality assessment of parameter estimates of signals embedded in Gaussian noise is extended in this paper over the case where useful signals are immersed in a mixture of compound Gaussian and Gaussian-distributed noises. The main problem here is that analytical expressions for distributions of such mixtures do not exist in most cases. Moreover, in some cases like K-distributed noise only, where closed-form expressions for the data distribution are available, the traditional Cramér-Rao bound does not exist. This makes the EL technique even more important for parameter estimation performance assessment. In this paper, for the so-called conditional model, we introduce test statistics whose distribution for the true (actual) parameters does not depend on these parameters and specifics of texture distribution, which makes them applicable for EL applications. We illustrate the utility of this EL technique by studying and predicting the performance breakdown of some direction of arrival estimators in a mixture of K-distributed and Gaussian noise.

I. INTRODUCTION

D IRECTION of arrival (DoA) estimation of multiple sources is a difficult problem that very rarely can be provided with the globally optimum solution, maximum likelihood (ML) in most cases. Indeed, only an exhaustive global search can guarantee strict optimality of its solution. For this reason, in most practical applications, different ML-proxy techniques are used. Most advanced are the so-called subspace-based techniques such as MUSIC [START_REF] Schmidt | Multiple emitter location and signal parameter estimation[END_REF] whereby instead of a global search over P DoA, one-dimensional (1D) search for the P maxima of the pseudo-spectrum function is required. In [START_REF] Stoica | MUSIC, maximum likelihood and Cramér-Rao bound[END_REF]- [START_REF] Stoica | MUSIC, maximum likelihood and Cramér-Rao bound: Further results and comparisons[END_REF] it was proven that these techniques may be treated as largesample equivalents of the accurate ML estimator for T → ∞

The associate editor coordinating the review of this manuscript and approving it for publication was Dr. Wei Liu. The work of O. Besson was supported by Délégation Générale de l'Armement/Mission pour la Recherche et l'Innovation Scientifique (DGA/MRIS) under Grant 2015.60.0090.00.470.75.01. Y. I. Abramovich is with the W R Systems, Ltd., Fairfax, VA 22030 USA (e-mail: yabramovich@wrsystems.com).

O. Besson is with the Institut Supérieur de l'Aéronautique et de l'Espace-Supaéro, Department of Electronics, Optronics, Signal, Toulouse 31055, France (e-mail: olivier.besson@isae-supaero.fr). B. A. Johnson is with the University of South Australia, Institute for Telecommunications Research, Mawson Lakes, SA 5085, Australia (e-mail: ben.a.johnson@ieee.org). independent and identically distributed (i.i.d.) Gaussian snapshots. In fact, it was proven in [START_REF] Stoica | MUSIC, maximum likelihood and Cramér-Rao bound[END_REF]- [START_REF] Stoica | MUSIC, maximum likelihood and Cramér-Rao bound: Further results and comparisons[END_REF] that asymptotic (T → ∞) DoA mean-square error (MSE) in this case tends to the Cramér-Rao bound. Similar results have been derived for other subspace techniques such as MODE [START_REF] Stoica | Novel eigenanalysis method for direction estimation[END_REF] or ESPRIT [START_REF] Roy | ESPRIT-Estimation of signal parameters via rotational invariance techniques[END_REF].

However, it was proven in [START_REF] Mestre | Modified subspace algorithms for DoA estimation with large arrays[END_REF] that, when M → ∞, T → ∞ and M/T → c, where M is the number of array elements, the localization function of MUSIC is not consistent and that MUSIC is not capable of consistently resolving sources within O(1/M ) angular distance. This inconsistency manifests itself in the socalled threshold area where, due to small sample support T and/or signal to noise ratio (SNR), the mean-square error of MUSIC is well above the CRB. This departure, which corresponds to the appearance of severely erroneous DoA estimates (outliers), occurs under threshold conditions (T, SNR) where the genuine ML estimator still provides CRB-consistent estimation accuracy [START_REF] Johnson | MUSIC, G-MUSIC and maximum likelihood performance breakdown[END_REF]. Mestre, based on random matrix theory, proposed in [START_REF] Mestre | Modified subspace algorithms for DoA estimation with large arrays[END_REF], [START_REF] Mestre | Improved estimation of eigenvalues and eigenvectors of covariance matrices using their sample estimates[END_REF], [START_REF] Mestre | On the asymptotic behavior of the sample estimates of eigenvalues and eigenvectors of covariance matrices[END_REF] an enhancement of MUSIC, namely G-MUSIC,which provides consistent estimates of the pseudospectrum function in the M → ∞, T → ∞, M/T → c regime. Johnson et al. [START_REF] Johnson | MUSIC, G-MUSIC and maximum likelihood performance breakdown[END_REF] investigated in detail the performance of G-MUSIC compared to that of ML, and provided explanations on why it improves over MUSIC without yet filling the gap between G-MUSIC-specific and ML-intrinsic threshold conditions.

For MUSIC, G-MUSIC as well as for other methods, much attention has been devoted on how to detect ("predict") the onset of the threshold regime, and how to possibly rectify ("cure") these method-specific outliers. A pivotal tool for such analysis of performance breakdown prediction is the so-called expected likelihood (EL) technique developed in [START_REF] Abramovich | Bounds on maximum likelihood ratio-Part I: Application to antenna array detection-estimation with perfect wavefront coherence[END_REF]- [START_REF] Abramovich | Expected likelihood support for deterministic maximum likelihood DOA estimation[END_REF]. The original EL approach was derived in the Gaussian case to test whether a candidate matrix, say R c , is a plausible covariance matrix for the observed data X ∈ C M ×T . It amounts to testing

H 0 : E XX H = R c against H 1 : E XX H = R
where R is unknown, and relies on the likelihood ratio for this test. A fundamental property is that LR(R 0 ), where R 0 is the true covariance matrix of X, has a distribution that depends only on M and T . Therefore, this distribution can be computed and used to assess the plausibility of R c . For DoA estimation problems where the covariance matrix depends on some parameter vector η (typically η includes DoA, sources power and white noise power), this approach enables one to assess the quality of an estimate η by computing LR(R(η)) and comparing it to LR(R(η 0 )) where η 0 is the true parameter vector. Naturally, the latter is not known a priori, but for the considered conditional [START_REF] Abramovich | Bounds on maximum likelihood ratio-Part I: Application to antenna array detection-estimation with perfect wavefront coherence[END_REF], [START_REF] Abramovich | Performance breakdown of subspacebased methods in arbitrary antenna arrays: GLRT-based prediction and cure[END_REF] and unconditional [START_REF] Abramovich | Expected likelihood support for deterministic maximum likelihood DOA estimation[END_REF] Gaussian models, the introduced LR(R(η 0 )) is shown to have a distribution that does not depend on η 0 , being specified by (M, T ) for unconditional models and by (M, T, P ) for conditional models. Specifically, it was demonstrated in [START_REF] Johnson | MUSIC, G-MUSIC and maximum likelihood performance breakdown[END_REF] that MUSIC-produced sets of DoA estimates that contain an outlier generate LR values that are statistically far smaller than the ones specified by the probability density function (p.d.f.) of LR(R(η 0 )). Thus the problem of statistically reliable prediction of MUSIC and G-MUSIC specific breakdown has been addressed in [START_REF] Johnson | MUSIC, G-MUSIC and maximum likelihood performance breakdown[END_REF]. Analysis of LR values generated by G-MUSIC demonstrated that G-MUSIC-specific outliers also produce LR values outside of the LR(R(η 0 )) distribution support. We also demonstrated that for certain small sample support T and/or SNR values, even genuine ML technique implemented via search for global likelihood function (LF) maximum starts to consistently produce estimation outliers. Yet, in this case, these ML-specific outliers are as likely as the true DoA values, which means the onset of the limits of the ML estimation paradigm, since the ML is no longer associated with the CRB accuracy in this ML-breakdown region. This EL analysis established the gap between the ML-proxy (MUSIC, G-MUSIC) specific and MLE-intrinsic threshold conditions and led to techniques capable of outliers rectification by driving the LR values produced by rectified DoA estimates into support of the LR(R(η 0 )) distribution. Majority of these results has been derived for Gaussian distributions although the EL approach was recently extended to complex elliptically symmetric distributions in [START_REF] Abramovich | Expected likelihood approach for covariance matrix estimation: Complex angular central Gaussian case[END_REF]- [START_REF] Besson | Regularized covariance matrix estimation in complex elliptically symmetric distributions using the expected likelihood approach-Part 2: The under-sampled case[END_REF], but mostly for the purpose of covariance matrix estimation.

In some applications, e.g., in HF direction finding applications, external noise is dominated by lightning strikes and being practically white for linear uniform arrays is strongly non-Gaussian. Also, in radar applications, the main source of external disturbance, in this instance clutter, is better described by a non-Gaussian distribution [START_REF] Conte | Characterisation of radar clutter as a spherically invariant process[END_REF]- [START_REF] Conte | Statistical tests for higher order analysis of radar clutter -Their analysis to L-band measured data[END_REF]. In [START_REF] Abramovich | Bounds for maximum likelihood regular and non-regular DoA estimation in K-distributed noise[END_REF], we addressed DoA estimation in pure K-distributed noise. The thermal Gaussian noise was neglected, on the ground that the non-Gaussian (external) to Gaussian (internal) noise ratio is generally high. Due to the availability of the p.d.f. of the observations in closedform, it was possible to derive the CRB (whenever it existed) as well as the ML estimator. However, an outcome of [START_REF] Abramovich | Bounds for maximum likelihood regular and non-regular DoA estimation in K-distributed noise[END_REF] was that internal noise should be taken into account even for high external to thermal noise ratio as, with very spiky clutter, some texture values can fall well below the thermal noise power. Therefore in [START_REF] Besson | Direction of arrival estimation in a mixture of K-distributed and Gaussian noise[END_REF] we addressed the more complicated problem of multiple DoA estimation of useful signals immersed into a mixture of K-distributed and Gaussian noises. A practical iterative DoA estimation scheme was introduced and compared to potential competitors such as the recently proposed robust G-MUSIC (RG-MUSIC) technique [START_REF] Couillet | Robust estimates of covariance matrices in the large dimensional regime[END_REF], [START_REF] Couillet | Robust spiked random matrices and a robust G-MUSIC estimator[END_REF]. In the large sample support or high SNR scenario, most methods were close to a lower bound derived under the (rather unrealistic) assumption of known texture of the compound Gaussian noise. The latter bound was achieved by the corresponding (clairvoyant) ML estimator which assumes all texture values known. In addition to being a not realistic assumption, in [START_REF] Besson | Direction of arrival estimation in a mixture of K-distributed and Gaussian noise[END_REF] the search for ML estimates was conducted in the close vicinity of the true DoA values, and therefore it does not reflect the threshold behavior of the genuine (global search) ML estimator. Moreover, as expected, for small sample support and/or SNR values, we observed strong departure from the clairvoyant lower bound. Yet, the behavior of all methods in this threshold area is not fully understood. Indeed, a major problem here is that the analytical distribution for such a noise mixture does not exist in closed-form, and therefore availability of a CRB or of the exact ML estimator as references is still an open problem in this case. In fact, in absence of an analytical expression for the likelihood function, the standard approach where brute-force global ML search for the optimal DoA estimates in the threshold area enables one to reveal the ML threshold conditions is not applicable. Therefore, an approach, similar to EL, is definitely required that would be capable of evaluating the proximity of the likelihood metric for the derived DoA estimates to the statistical values of this metric calculated at the true DoA values. This is the aim of this paper.

Specifically, the paper is organized as follows. In Section II, we present the model at hand, briefly review the DoA estimators derived in [START_REF] Besson | Direction of arrival estimation in a mixture of K-distributed and Gaussian noise[END_REF] and we propose a new estimator which enables one to relax some assumptions made in [START_REF] Besson | Direction of arrival estimation in a mixture of K-distributed and Gaussian noise[END_REF]. Then, we introduce two test statistics which meet the requirements for EL application, i.e., their distributions, when evaluated at the true parameters, depend only on (M, T, P ). Therefore, they allow for assessment of the estimators introduced above, especially in the threshold area. This is the aim of the simulations of Section III where we evaluate the mean-square error of the various estimators and show that outliers can be safely detected by the EL approach. The analysis enables one to reveal the estimators-specific as well as the ML-intrinsic breakdown conditions.

II. DATA MODEL, DOA ESTIMATORS AND EXPECTED LIKELIHOOD APPROACH

A. Data Model

The data model is essentially that of [START_REF] Besson | Direction of arrival estimation in a mixture of K-distributed and Gaussian noise[END_REF] and is given by

x t = A(θ)s t + √ τ t n t + σ w w t ( 1
)
where

θ = θ 1 • • • θ P T is the vector of the DoA and A(θ) = a(θ 1 ) • • • a(θ P )
is the manifold matrix of the M -element array. We assume a conditional model for which the emitted waveforms s t are treated as deterministic unknowns. n t and w t are independent and identically distributed vectors drawn from n t ∼ CN (0, I) and w t ∼ CN (0, I). σ 2 w stands for the thermal noise power. τ t is a positive variable which we choose to treat as deterministic unknown, similarly to S. The choice of a conditional model rather than an unconditional model where S and τ would be considered random with a prior distribution is due to the fact that these prior distributions may not be completely known and that derivation of maximum likelihood estimators is simpler in the conditional model.

B. DoA Estimation

For the sake of clarity, we first provide a quick overview of the two main estimators derived in [START_REF] Besson | Direction of arrival estimation in a mixture of K-distributed and Gaussian noise[END_REF] and then we propose a new exact ML estimator which relaxes some assumptions made in [START_REF] Besson | Direction of arrival estimation in a mixture of K-distributed and Gaussian noise[END_REF]. Under the stated assumptions, the p.d.f. of the data matrix

X = x 1 • • • x T is given by p(X; θ, τ , S) ∝ T t=1 τ t + σ 2 w -M exp -τ t + σ 2 w -1 x t -A(θ)s t 2 (2)
where

∝ means proportional to, τ = τ 1 • • • τ T T and S = s 1 • • • s T .
In [START_REF] Besson | Direction of arrival estimation in a mixture of K-distributed and Gaussian noise[END_REF] we first investigated a clairvoyant maximum likelihood estimator of θ, assuming that τ and σ 2 w are known. This clairvoyant ML estimator of θ was obtained as

θ ML |τ ,σ 2 w = arg min θ T t=1 x H t P ⊥ A(θ) x t τ t + σ 2 w ( 3
)
where H . In [START_REF] Besson | Direction of arrival estimation in a mixture of K-distributed and Gaussian noise[END_REF],

P ⊥ A(θ) = I -A(θ) A(θ) H A(θ) -1 A(θ)
only a local search for the maximum of the concentrated likelihood function in (3) was implemented, meaning that θ ML |τ ,σ 2 w was sought in close vicinity of the true DoA θ 0 . A lower bound, namely the averaged-over-τ CRB conditioned on τ , i.e., E {CRB(θ|τ )} was also derived. With the local search of the maximum implemented, it was observed that the meansquare error (MSE) was very close to this bound, whatever SNR. However, the performance of the genuine ML estimator, which consists of a global search for the maximum of T t=1 (τ t + σ 2 w ) -1 x H t P ⊥ A(θ) x t was not studied in [START_REF] Besson | Direction of arrival estimation in a mixture of K-distributed and Gaussian noise[END_REF], which would have revealed the ML-intrinsic performance breakdown conditions. We will fill this gap in the next section.

In order to relax the assumption of known τ t , we also considered an approximate maximum likelihood estimation of both the unstructured steering matrix A and τ , assuming only σ 2 w is known. This method relies on the fact that, for a given A, the MLE of τ t is

τ ML t|A,σ 2 w = max M -1 x H t P ⊥ A x t , σ 2 w -σ 2 w (4)
while, for known τ , the MLE of A is

A ML |τ ,σ 2 w = P P T t=1 x t x H t τ t + σ 2 w ( 5
)
where P P (.) stands for the P principal subspace of the matrix between parentheses. The approximate maximum likelihood (AML) estimator of A and τ was thus implemented through the iterative procedure described in Table I. We let A AML |σ 2 w denote the value obtained at the end of the simulations. The DoA were then estimated as

θ AML |σ 2 w = MUSIC A AML |σ 2 w (6)
where MUSIC [.] stands for the conventional MUSIC algorithm applied to a M × P matrix whose columns form a basis for the signal subspace.

In the present paper, we go one step further by assuming that neither τ nor σ 2 w are known and we derive the maximum likelihood estimator of θ and γ t = τ t + σ 2 w . Indeed, when σ 2 w is unknown, one should consider γ t = τ t + σ 2 w as a whole unknown since the p.d.f. of the data matrix X is given by p(X; θ, γ, S)

∝ T t=1 γ -M t exp -γ -1 t x t -A(θ)s t 2 (7) where γ = γ 1 • • • γ T T . It is readily verified that max S p(X|A(θ), γ, S) ∝ T t=1 γ -M t exp -γ -1 t x H t P ⊥ A(θ) x t . (8)
Moreover,

max γ t γ -M t exp -γ -1 t x H t P ⊥ A(θ) x t ∝ x H t P ⊥ A(θ) x t -M (9)
so that the ML estimator of θ is obtained as

θ ML = arg max θ T t=1 x H t P ⊥ A(θ) x t -M = arg max θ - T t=1 log x H t P ⊥ A(θ) x t . ( 10
)
In the simulations section, we will compare the performance of θ ML to that of θ ML |τ ,σ 2 w .

C. Expected Likelihood Approach

In order to assess the various estimators above and others, i.e., in order to detect potential outliers, we suggest to resort to the EL principle and, consequently, we look for a test statistic (actually a likelihood ratio) whose distribution, when evaluated at the true parameters, does not depend on S or τ or θ 0 but only on known parameters. Then, to assess the quality of any given estimate, the test statistic evaluated at the estimated parameters will be checked against the test statistic at the true parameters. If a discrepancy is observed, typically when the LR for the estimated parameters fall outside the support of the LR for the true parameters, then the estimates are deemed erroneous. Towards this end, one could think of using the likelihood function in [START_REF] Mestre | On the asymptotic behavior of the sample estimates of eigenvalues and eigenvectors of covariance matrices[END_REF], evaluated at some estimate θ, for assessment of the latter, by comparison with the LF evaluated at the true DoA θ 0 . However, when evaluated at the true steering matrix A 0 , one has

x H t P ⊥ A 0 x t = [ √ τ t n t + σ w w t ] H P ⊥ A 0 [ √ τ t n t + σ w w t ] = [ √ τ t n t + σ w w t ] H U 0 U H 0 [ √ τ t n t + σ w w t ] d = (τ t + σ 2 w ) wH t U 0 U H 0 wt (11)
where d = means "distributed as". In the previous equation, we used the fact that, when τ t is deterministic, √ τ t n t + σ w w t ∼ CN 0, (τ t + σ 2 w )I M and hence it is distributed as τ t + σ 2 w wt where wt ∼ CN (0, I M ). Moreover U 0 is an M × (M -P ) semi-unitary matrix whose columns form an orthonormal basis for R(A 0 ) ⊥ , the subspace orthogonal to A 0 , and hence U H 0 wt ∼ CN (0, I M -P ). This likelihood function depends on τ t , and hence its distribution is dependent on the statistics of τ t . To get rid of this dependency, one might consider, for any estimate θ and corresponding basis U( θ) for the associated noise subspace, the following data

z t ( θ) = U( θ) H x t U( θ) H x t (12)
which, for the true θ 0 , follows a complex angular central Gaussian (ACG) distribution with parameter matrix I M -P , since

z t (θ 0 ) = U H 0 x t U H 0 x t d = U H 0 wt U H 0 wt
.

In other words, z t (θ 0 ) is uniformly distributed on the unit sphere. It follows that quality of an estimate θ may be assessed by proximity of the scatter matrix of z t ( θ) to the identity matrix. This suggests the use of the traditional sphericity test [START_REF] Muirhead | Aspects of Multivariate Statistical Theory[END_REF] ST θ = det(R( θ))

(M -P ) -1 Tr{R( θ)} M -P (13)
where

R( θ) = T -1 T t=1 z t ( θ)z H t ( θ). ( 14
)
Note that, since z t ( θ) is unit norm, ST ( θ) ∝ det(R( θ)).

The sphericity test measures the spreading of the eigenvalues of the matrix R( θ), whose sum is equal to one. The maximum is achieved when all eigenvalues are equal, which amounts to R( θ) being the identity matrix. The suggested here "expected likelihood" test is a monotonic function of det(R( θ)) where the threshold for acceptance or rejection of θ is defined by the distribution of det(R(θ 0 )), which is that of det(T -1 T t=1 z t (θ 0 )z H t (θ 0 )) with z t (θ 0 ) uniformly distributed on the complex sphere. Note that, in the absence of an analytical expression for the p.d.f. of this statistic, the threshold is computed from Monte-Carlo simulations. The introduced test checks the "whiteness" of data after projection onto the subspace orthogonal to A( θ). The EL sphericity test meets all requirements: its distribution, when evaluated at the true DoA θ 0 , does not depend on the DoA nor on the statistics of the texture distribution, and is only specified by M , P and T .

The sphericity test, which is meant at testing the spreading of the eigenvalues of R( θ) can also be viewed a generalized likelihood ratio for testing whether the covariance matrix of a Gaussian distributed data is the identity matrix. Since, z t (θ 0 ) follows a complex angular central Gaussian (ACG) distribution with parameter matrix I M -P , this suggests using the likelihood ratio for testing whether the scatter matrix of z t ( θ) is the identity matrix, considering that z t ( θ) follows an ACG distribution. Such generalized likelihood ratio test for ACG distributions was indeed derived in [START_REF] Abramovich | Regularized covariance matrix estimation in complex elliptically symmetric distributions using the expected likelihood approach-Part 1: The oversampled case[END_REF], [START_REF] Besson | Regularized covariance matrix estimation in complex elliptically symmetric distributions using the expected likelihood approach-Part 2: The under-sampled case[END_REF]. Let Z = z 1 • • • z T be a (M -P ) × T matrix where each z t follows a complex ACG distribution with scatter matrix Σ 0 . Let Σ c be a candidate for this scatter matrix. Then the LR for testing Σ c is

LR ACG = det(Σ ML Σ -1 c ) T T t=1 z H t Σ -1 c z t z H t Σ -1 ML z t -(M -P ) (15)
where Σ ML obeys

Σ ML = M -P T T t=1 z t z H t z H t Σ -1 ML z t . ( 16
)
As shown in [START_REF] Abramovich | Regularized covariance matrix estimation in complex elliptically symmetric distributions using the expected likelihood approach-Part 1: The oversampled case[END_REF], [START_REF] Besson | Regularized covariance matrix estimation in complex elliptically symmetric distributions using the expected likelihood approach-Part 2: The under-sampled case[END_REF], this LR, when evaluated at the true scatter matrix Σ 0 (equal to I M -P in our case) has a distribution that depends on M , P and T only. Therefore, such a likelihood ratio, applied to z t ( θ), can serve as a means to detect if the data is white after projection. The difference compared to the sphericity test is that these LR are not based on the same assumptions. Note also that LR ACG is more complicated to compute than the sphericity test.

III. SIMULATIONS

We consider a uniform linear array of M = 20 elements spaced a half-wavelength apart and a scenario with two sources, whose DoA are θ 1 = 16 • and θ 2 = 18 • . The waveforms are generated from a Gaussian distribution with covariance matrix P s I 2 . The compound-Gaussian noise follows a Weibull distribution: the τ t are Gamma distributed with shape parameter ν and scale parameter ν -1 so that E {τ t } = 1. We set ν = 0.2 in the sequel. The total noise power is thus 1 + σ 2 w and the non Gaussian to Gaussian noise ratio is simply nGGNR = 10 log 10 σ -2 w . The signal to noise ratio is defined as SNR = 10 log 10 (P s /(1 + σ 2 w )). 10 5 Monte-Carlo simulations were run to estimate the meansquare error (MSE) of the various estimators.

As for θ ML |τ ,σ 2 w and θ ML , the search for the maximum of the corresponding likelihood function was carried out using two protocols:

1) In a first implementation, the region {-π/2 ≤ θ 1 < π/2; θ 1 < θ 2 )} was sampled on a grid with step Δθ = 1 • and the search for the maximum of the LF was carried out on this grid, yielding an initial estimate (θ

2 ). Then, the search was refined on [θ (0)

1 -Δθ/2, θ (0) 1 + Δθ/2] × [θ (0) 2 -Δθ/2, θ (0) 2 + Δθ/2]
. This method is referred to as "global" in the figures, since it provides the global maximum of the LF. 2) A second implementation consists in searching the maximum of the LF in the vicinity of the true DoA: more precisely, the maximum was constrained to be at most a half beam-width apart from θ 0 . This method is referred to as "local" in the figures. We also consider RG-MUSIC [START_REF] Couillet | Robust estimates of covariance matrices in the large dimensional regime[END_REF], [START_REF] Couillet | Robust spiked random matrices and a robust G-MUSIC estimator[END_REF] which computes a consistent estimate (in the random matrix theory sense) of the noise projection matrix from the eigenvalue decomposition of

RRG-MUSIC = lim

k →∞ Rk with Rk+1 = 1 T T t=1 u 1 M x H t R-1 k x t x t x H t .
In the simulations u(x) = (1 + α)/(α + x) with α = 0.01.

A. Mean-Square Error of DoA Estimates

We first focus on evaluating the mean-square error of the various estimators, with two main objectives: 1) reveal the MLE-intrinsic performance breakdown conditions i.e., the SNR conditions under which the global search for the maximum of the likelihood function completely fails in providing reliable estimates and 2) compare θ ML |τ ,σ 2 w of [START_REF] Besson | Direction of arrival estimation in a mixture of K-distributed and Gaussian noise[END_REF] which assumes both τ t and σ 2 w are known to the present θ ML which does not make this assumption. Fig. 1 displays the MSE versus SNR. Two main observations can be made. First, θ ML |τ ,σ 2 w and θ ML provide the same MSE, whatever the type of implementation, local or global search for the maximum. This means that θ ML is an interesting option as it does not require as many assumptions as θ ML |τ ,σ 2 w . Secondly, a significant difference is seen between the two implementations. When a local search is performed, both θ ML |τ ,σ 2 w and θ ML stick to E {CRB(θ|τ )} whatever SNR, which confirms that the latter is a local bound. In contrast, with the global search for the maximum, under some SNR threshold, both θ ML |τ ,σ 2 w and θ ML significantly depart form the bound revealing the intrinsic-ML performance breakdown conditions. As was already observed for MUSIC and G-MUSIC, RG-MUSIC and θ AML |σ 2 w begin to produce bad results at SNR where the global MLE is still consistent, and the difference can be significant, e.g., 20 dB at T = 40 and nGGNR = 30 dB. As will be shown below, the distribution of the sphericity test or the likelihood ratio enables one to capture this phenomenon.

B. Distribution of the Likelihood Ratios and Detection of Outliers

In order to understand these different breakdown conditions, the sphericity test (ST) and LLR ACG = log LR ACG are instrumental, as illustrated in Figs. 2 and3 where we plot their distributions, for three values of SNR: SNR = -18 dB where even ML breaks, SNR = -6 dB where ML produces consistent estimates while RG-MUSIC and θ AML |σ 2 w do not, and SNR = 12 dB where all algorithms achieve the lower bound. It is clear that, when even MLE breaks down, it produces estimates whose ST or LLR is commensurate with that at the true DoA, which is a symptom of complete failure. On the other hand, for SNR = -6 dB, it is clearly seen that only RG-MUSIC and θ AML |σ 2 w produce DoA estimates whose ST or LLR is outside the support of the p.d.f. for θ 0 . Finally, at SNR = 12 dB, all estimates produce ST or LLR commensurate with those at θ 0 .

We now illustrate how the EL principle, through examination of ST and LLR ACG enable one to detect outliers, i.e., severely erroneous DoA estimates. Towards this end, we proceeded in the following way. For each method, among the 10 5 Monte-Carlo simulations, we decided that an estimate θ is a "true outlier" if θθ 0 2 is superior to 100 times the corresponding Cramér-Rao bound. Next, the sphericity test or the log likelihood ratio were computed for this θ: whenever they are below a threshold, θ is declared an outlier according to the EL principle. The probability of correct outlier detection is the ratio of the number of declared outliers in the set of true outliers to the number of true outliers. The thresholds of the ST or the LR were set to ensure a probability of false alarm equal to 10 -2 when ST or LLR are calculated with the true scatter matrix I M -P . Figs. 4 and5 plot the actual percentage of outliers and the probability of detecting them by comparing the ST or LLR to a threshold. As can be observed, for very low SNR where even ML breaks down, there is no way to detect these outliers. However, in the threshold area, where ML is accurate but RG-MUSIC or θ AML |σ 2 w produce erroneous estimates, the latter can be rather reliably detected using the EL principle. We notice that LLR ACG is slightly better than the sphericity test. However, the difference is not important and ST is easier to obtain than LLR ACG .

To illustrate further how this comparison of ST or LLR to a threshold can help, we plot in Figs. 6 and7 the MSE obtained when all estimates below the LLR threshold have been removed. Clearly, an improvement is observed, especially in the threshold area where the MSE is noticeably decreased. However, the gap between θ AML |σ 2 w or RG-MUSIC and θ ML is not completely filled with the selected false alarm threshold. By raising this threshold, we may improve MSE of the selected estimates further at expense of discarding more legitimate estimates.

IV. CONCLUSIONS

In this paper we continued our investigation into the DoA estimation problem for signals immersed in a mixture of external compound Gaussian noise and internal Gaussian noise, with the main focus on studying the threshold performance of DoA estimation routines. Since the p.d.f. of the noise mixture is not available, we resorted to the maximum likelihood methodology by treating random texture values as deterministic unknown parameters. The conditional likelihood function and associated likelihood ratio which tests sphericity of the projected out snapshots have been introduced for DoA estimates quality assessment using the expected likelihood principle. Specifically, the introduced sphericity test and likelihood ratio for the true DoA is described by distribution that does not depend on DoA nor on texture distribution, but is fully specified by antenna dimension, number of sources and sample support. Efficiency of estimation-method-specific outliers detection and ML-intrinsic breakdown prediction by the EL approach has been demonstrated.

Fig. 1 .

 1 Fig. 1. Mean square error of estimators versus SNR. ν = 0.2.

Fig. 2 .

 2 Fig. 2. Probability density function of sphericity test and LLR ACG for various SNR. T = 40, ν = 0.2 and nGGNR = 30 dB.

Fig. 3 .

 3 Fig. 3. Probability density function of sphericity test and LLR ACG for various SNR. T = 80, ν = 0.2 and nGGNR = 30 dB.

Fig. 4 . 2 w.

 42 Fig. 4. Probability of detecting outliers using the EL approach for RG-MUSIC and θ AML |σ 2 w . P f a = 10 -2 , T = 40, ν = 0.2 and nGGNR = 30 dB.

Fig. 5 . 2 w.

 52 Fig. 5. Probability of detecting outliers using the EL approach for RG-MUSIC and θ AML |σ 2 w . P f a = 10 -2 , T = 80, ν = 0.2 and nGGNR = 30 dB.

Fig. 6 .

 6 Fig. 6. MSE of RG-MUSIC and θ AML |σ 2 w with/without thresholding . P f a = 10 -2 , T = 40, ν = 0.2 and nGGNR = 30 dB.

Fig. 7 . 2 w

 72 Fig. 7. MSE of RG-MUSIC and θ AML |σ 2 w with/without thresholding . P f a = 10 -2 , T = 40, ν = 0.2 and nGGNR = 60 dB.

TABLE I APPROXIMATE

 I MAXIMUM LIKELIHOOD ESTIMATION OF A AND τ

FOR KNOWN σ 2 w


Mohamed L Derouiche
email: mohamedlotfi.derouiche@enit.rnu.tn

Soufiene Bouall Ègue
email: soufiene.bouallegue@issig.rnu.tn

Joseph Hagg Ège
email: joseph.haggege@enit.rnu.tn

Guillaume Sandou
email: guillaume.sandou@centralesupelec.fr

Rapid Model Predictive Control Prototyping with LabVIEW/CDSim and CompactRIO Target

Keywords: model predictive control, computer aided design, PIL co-simulation, rapid prototyping, LabVIEW/CDSim, NI Compact RIO-9082 platform, MAGLEV system

In this paper, an advanced Computer Aided Design (CAD) methodology for the Process-In-the-Loop (PIL) cosimulation and rapid prototyping of model predictive controllers is proposed and successfully implemented using the NI CompactRIO-9082 RT target and a host PC. The developed software (SW) and hardware (HW) co-design platform is based on the Control Design and Simulation (CDSim) module of LabVIEW environment and the Network Streams data communication protocol. The designed LabVIEW-based MPC algorithm as well as the dynamic model of the controlled plant are implemented as VI under the cRIO-9082 target and the host PC, respectively. This hardware model will be deployed on the CompactRIO Real-Time (RT) target within a PIL co-simulation framework. The proposed NI CompactRIO-9082 based CAD approach for prototyping and implementation of MPC algorithms is applied to the position control of a Magnetic Levitation system (MAGLEV). All obtained SW/HW co-simulation results are compared and discussed in order to improve the effectiveness of the proposed MPC co-design methodology.

I. INTRODUCTION

Model Predictive Control (MPC) is an advanced design approach that focuses on constructing and optimizing feedback controllers at each discrete-time instant [START_REF] Rawlings | Model Predictive Control: Theory and Design[END_REF]- [START_REF] Wang | Model Predictive Control System Design and Implementation Using MATLAB[END_REF]. Such predictive controllers can adjust the control action before a change in the output setpoint actually occurs. The ability of handling operational constraints in an explicit manner is one of the main reasons for the popularity and successful of MPC approaches in industrial applications especially for electrical motors drives, power systems, chemical process and so on [START_REF] Mayne | Model predictive control: Recent developments and future promise[END_REF], [START_REF] Joe Qin | A survey of industrial model predictive control technology[END_REF].

The implementation of model predictive controllers for complex dynamic systems requires usually non negligible software processing power, resources and hardware materials. The complexity of MPC algorithms (on-line constrained optimization, prediction, estimation, etc.) makes difficult the practical implementation of these digital control laws. Such an implementation becomes increasingly complex with the complexity of the controlled dynamic plants. The textual programming based method of MPC algorithms implementation, such as with C, C++ and VHDL languages, is a tedious and error-prone task [START_REF] Bouafassa | Design and real time implementation of single phase boost power factor correction converter[END_REF]- [START_REF] Bolognani | Design and Implementation of Model Predictive Control for Electrical Motor Drives[END_REF]. This may present a real brake to the design of increasingly complex MPC controllers. However, for rapid prototyping of these predictive control structures purposes, the use of high level development tools such as graphical code generators such as the LabVIEW Control Design and Simulation module (LabVIEW/CDSim) [START_REF]National Instruments, LabVIEW Control Design User Manual[END_REF] and the hardware Compact Reconfigurable Input/Output (cRIO) platforms [START_REF]NI LabVIEW for Compact RIO Developer's Guide[END_REF], [START_REF]CompactRIO NI cRIO-9081/9082[END_REF] is a promising solution for the easy and fast MPC algorithms implementation. The Control Design and Simulation software of LabVIEW, associated to a NI multi-core cRIO RT target and a host PC, presents an efficient tool for Processor-In-the-Loop (PIL) and Hardware-In-the-Loop (HIL) co-simulations and rapid prototyping [START_REF] Boumazbar | Co-simulation and Rapid Prototyping of Fuzzy Supervised PID controllers based on FPGA[END_REF], [START_REF] Boumazbar | Co-simulation and Rapid Prototyping of Fuzzy Controllers based on FPGA Nexys2 Board[END_REF]. These software/hardware (SW/HW) environments and resources will be used together to develop a powerful PIL co-design platform for MPC implementation framework.

The remainder of this paper is organized as follows. In Section II, the model predictive control concept is introduced. A LabVIEW/CDSim-based workflow to design and implement MPC controllers is particularly described. Section III presents the proposed NI cRIO-9082 based CAD methodology for PIL co-simulation and rapid prototyping of MPC controllers. An application to rapid control prototyping of a magnetic levitation system is investigated in Section IV. All SW/HW cosimulation results, obtained in the proposed CAD methodology with LabVIEW/CDSim and cRIO tools, are presented and discussed. Section V concludes this paper.

II. LABVIEW/CDSIM-BASED MODEL PREDICTIVE CONTROL

This section provides information about using the Lab-VIEW Control Design and Simulation Module to design and implement MPC controllers. Basic concepts about MPC design are firstly presented.

A. Basic concepts of MPC design

As mentioned by J.A. Rossiter in [START_REF] Rossiter | Model Based Predictive Control: A practical Approach, Control Series[END_REF], MPC approach reflects human behavior whereby we select control actions which we think will lead to the best predicted outcome over some limited horizon. To make this selection, an internal model of the process in question is used. Decisions are constantly updated as new observations become available. So, an MPC algorithm consists mainly of a cost function to be minimized, operational constraints and model of the controlled plant [START_REF] Rawlings | Model Predictive Control: Theory and Design[END_REF]- [START_REF] Wang | Model Predictive Control System Design and Implementation Using MATLAB[END_REF]. The prediction horizon N p , used in MPC approach, is the number of samples in the future during which the feedback controller predicts the plant output. The control horizon N c is the number of samples within the prediction horizon during which the predictive controller can affect the control action as depicted in Fig. 1. The value specified for this horizon must be less than the one used for the prediction horizon. So, the main elements of the discrete-time model predictive control are depicted in Fig. 2. The plant input, the controlled output and the reference trajectory are denoted by u ∈ R, y ∈ R, and r ∈ R, respectively. The plant model is employed to calculate output predictions up to samples of prediction horizon in the future. The optimization algorithm is aimed at determining the control sequence given by {u (k -1 + i) , i = 1, 2, . . . , N c }. This sequence minimizes the specified cost function [START_REF] Rawlings | Model Predictive Control: Theory and Design[END_REF] in the MPC design formalism, subject to problem constraints on the input, change in input and output of the plant. The optimization algorithm assumes that u (k - For SISO systems, and when omitted the term related on the control action error, the optimization problem, which penalizes tracking errors and control variations and often used in the MPC framework, is given as follows [START_REF] Rossiter | Model Based Predictive Control: A practical Approach, Control Series[END_REF]:

1 + i) = u (k + N c -1) for N c < i ≤ N p .
min ∆U J = Np i=1 [r (k + i) -ŷ (k + i|k)] 2 +λ Nc i=1 [∆u (k -1 + i)] 2 (1
) where ∆U = [∆u (k) , ∆u (k + 1) , . . . , ∆u (k + N c -1)] T is the vector of decision variables and ∆u (k) = u (k) -u (k -1). The design parameter λ > 0, called weighting factor, is adjusted to achieve a compromise between minimizing the output tracking error and minimizing variations on the control signal. In addition, ŷ (k + i|k) denotes the prediction of the output at instant k + i on the basis of the measured state x (k) ∈ R.

From an implementation point of view, only the first element of the optimized control sequence is applied to the plant and the control input is updated at each sampling time. The optimization process is repeated at the next sampling instant, on the basis of the measured state vector x (k). In addition to the weight coefficient λ > 0 (or weighting matrices in the case of MIMO systems) in the cost function, operational constraints on the control action, its rate of change and plant output signals limitations can be defined, respectively, as given by Equations ( 2), ( 3) and ( 4):

u min ≤ u (k) ≤ u max ( 2
)
∆u min ≤ ∆u (k) ≤ ∆u max (3)
y min ≤ y (k) ≤ y max (4)

B. MPC design with LabVIEW/CDSim

To create an MPC controller under LabVIEW/CDSim environment, the CD Create MPC Controller VI of Fig. 3 is used. This VI bases the MPC controller on a state-space model of the controlled plant. The CD Create MPC Controller VI is also used to specify constraints for the predictive controller design. MPC constraints are specified by either the dual optimization or the barrier function method. These handling constraints methods are manually selected to use by the polymorphic instances of such a VI. More details about these methods are given in [START_REF]National Instruments, LabVIEW Control Design User Manual[END_REF]. The constraints for MPC design problem are specified in the MPC Constraints parameter input of the CD Create MPC Controller VI. All these constraints can be updated in run time [START_REF]National Instruments, LabVIEW Control Design User Manual[END_REF].

Barrier type approach of MPC handling constraints specifies any constraints on the quadratic programming optimization algorithm using the barrier function method. Each parameter specifies a minimum limit, maximum limit, tolerance, or penalty on the initial or final control action u (k), plant output y (k), or rate of change in control action ∆u (k). When specifying only a minimum or a maximum limit for a parameter, LabVIEW assumes a constant profile for that parameter. Dual MPC constraints method specifies any constraints on the quadratic programming optimization algorithm using the dual optimization method.

Finally, the created MPC controller is now implemented either in a simulation or a real-world scenario while using the CD Implement MPC Controller VI of Fig. 4 within a timed or simulation loops [START_REF]National Instruments, LabVIEW Control Design User Manual[END_REF]. This VI needs, in addition to the created MPC controller, information about profiles of the input setpoints (Control Action Reference Window) and output setpoints (Output Reference Window), disturbances (Disturbance Window) and mainly the measured output y (k) of the controlled plant, to calculate the control action u (k) to apply to the plant along the control horizon at time k. According to Fig. 4, the CD Implement MPC Controller VI returns mainly the information related to the control action necessary to react to the change in the output setpoint profile, the predicted output of the plant along the prediction horizon and the rate of change in control action. In all cases, it is possible to provide setpoint and disturbance profile information either in advance of controller execution or dynamically during this execution [START_REF]National Instruments, LabVIEW Control Design User Manual[END_REF].

III. PROPOSED SW/HW CO-SIMULATION APPROACH FOR MPC PROTOTYPING

A. CAD methodology for rapid MPC prototyping

A powerful PIL co-simulation platform for MPC algorithms verification and prototyping must make easy the practical implementation of such control laws with the same hardware target and software tools [START_REF] Boumazbar | Co-simulation and Rapid Prototyping of Fuzzy Supervised PID controllers based on FPGA[END_REF], [START_REF] Boumazbar | Co-simulation and Rapid Prototyping of Fuzzy Controllers based on FPGA Nexys2 Board[END_REF]. For this purpose, the principle of our proposed CAD methodology for rapid The developed SW/HW solution is based mainly on the use of a multi-core CompactRIO-9082 platform from National Instruments [START_REF]NI LabVIEW for Compact RIO Developer's Guide[END_REF], [START_REF]CompactRIO NI cRIO-9081/9082[END_REF]. Well-suited to complex processing and real-time computing, this RT target is associated to a host PC equipped with LabVIEW/CDSim environment. This advanced control platform is connected to the host PC via an Ethernet connection during prototyping phase. The cRIO-9082 device operates autonomously to execute a control LabVIEW project deployed on its RT dual-core processor.

Thanks to its powerful tools, LabVIEW software simplifies construction and prototyping of designed control systems and provides the ability to implement a variety of control algorithms, in particular those of MPC strategy. The integrated Control Design & Simulation Module of LabVIEW presents a powerful tool to implement, prototype and test such algorithms [START_REF]National Instruments, LabVIEW Control Design User Manual[END_REF].

B. NI Compact RIO-9082 based implementation

The main component of the proposed PIL platform for rapid prototyping of MPC design is the NI Compact RIO-9082 target, Fig. 6. This high-performance multi-core NI cRIO-9082 system provides advanced Intel Core i7 dual-core processor with 1.33 GHz frequency speed, 32 GB nonvolatile storage and 2 GB DDR3 800 MHz RAM [START_REF]NI LabVIEW for Compact RIO Developer's Guide[END_REF], [START_REF]CompactRIO NI cRIO-9081/9082[END_REF]. This hardware target provides also a built-in VGA display output for an integrated user interface and the option to use a Microsoft Windows Embedded Standard 7 (WES7) or LabVIEW Real-Time Operating System (OS).

The NI cRIO-9082 system is built around an 8-slot Spartan-6 LX150 FPGA chassis for custom I/O timing, control, and processing. This chassis contains 1 MXI-Express, 4 USB Hi-Speed, 2 Gigabit Ethernet, and 2 serial ports for connectivity and expansion. An operating temperature range of 0 to 55 • C is tolerated for this multi-core RT target. The NI cRIO-9082 device offers the widest array of connectivity and expansion options available in the NI Com-pactRIO platform, including the high-bandwidth and lowlatency MXI-Express bus for expansion using the multi-slot MXI-Express RIO chassis. The increased processing power of the NI cRIO-9082 target makes it well suited to perform the advanced processing tasks required by complex and hard applications such as rapid control prototyping and HIL cosimulation.

The LabVIEW Real-Time workflow is chosen to take advantage of deterministic execution and the highest degree of reliability in continuous operation environments [START_REF]CompactRIO NI cRIO-9081/9082[END_REF]. However, the WES7 based one takes advantage of the extensive Windows ecosystem of software and display capabilities made possible by LabVIEW software.

IV. APPLICATION TO MPC PROTOTYPING OF A MAGLEV SYSTEM

A. Plant description and modeling

The magnetic levitation system (MAGLEV) of Feedback company [START_REF]-942S[END_REF], which is depicted in Fig. 7, is used as a process example to implement the proposed PIL co-simulation approach.

In required range of operation, the distance h of the suspended sphere is given by the infrared photo-sensor voltage y as follows [START_REF]-942S[END_REF]:

y = γh + y 0 ( 5
)
where γ is a positive gain depending on the position sensor, and y 0 is the offset voltage such that y ∈ [-2V, +2V ].

The coil current is regulated by an inner control loop within the driver block. Its characteristic is linearly related to the input voltage u as follows, neglecting its high frequency dynamics: where ρ > 0 is the coil resistor and i 0 > 0 is the offset value of current.

i = ρu + i 0 (6)
The working excursion of u is limited between -3V, corresponding to a null coil current, and +5V that defines the saturation value [START_REF]-942S[END_REF]- [START_REF] Fama | Predictive Control of a Magnetic Levitation System with Explicit Treatment of Operational Constraints[END_REF]. For this process, the predictive control of the suspended sphere position is investigated. The dynamics of the vertical movement of such a sphere is modeled as follows [START_REF]-942S[END_REF], [START_REF] Santos | Robust Model Predictive Control for a Magnetic Levitation System Employing Linear Matrix Inequalities[END_REF], [START_REF] Fama | Predictive Control of a Magnetic Levitation System with Explicit Treatment of Operational Constraints[END_REF]:

m d 2 h dt 2 = mg -K i 2 h 2 (7)
where K is an electromechanical conversion gain depending on the MAGLEV system, m is the mass of the sphere, g is the acceleration of gravity, and i is the coil current.

According to the given sensor and current driver characteristics ( 5) and ( 6), the equation ( 7) can be re-written as:

m d 2 y dt 2 = γmg -K (ρu + i 0 ) 2 γ 3 (y -y 0 ) 2 (8)
Taking x = [ y ẏ ] T as state vector, the following state space representation of the studied system is obtained:

ẋ1 = x 2 ẋ2 = γg -K (ρu+i0) 2 γ 3 m(x1-y0) 2 (9)
Finally, while using the values for the physical model parameters, a linear discrete-time model of the MAGLEV plant, where a zero-order-hold with 5 ms sampling period was adopted at the input of the system [START_REF] Santos | Robust Model Predictive Control for a Magnetic Levitation System Employing Linear Matrix Inequalities[END_REF], [START_REF] Fama | Predictive Control of a Magnetic Levitation System with Explicit Treatment of Operational Constraints[END_REF], is given by the following state-space representation:

   x k+1 = 1.0108 0.0050 4.3185 1.0108

x k + -0.0142 -5.6779 u k y k = [ 1 0 ] x k (10)

B. Implementation and practical results

In this section, the proposed PIL co-simulation approach for rapid prototyping of the MPC of MAGLEV system is described. All implementation and practical results will be given and discussed.

To implement the proposed PIL co-simulation of MPC design, according to the principle Fig. 5, two separate VIs are developed. The first one, dedicated to the cRIO-9082 target deployement, is used to configure the MPC controller by introducing the prediction model, input and output constraints and weight parameters of the cost function (see Fig. 8). This VI Fig. 8: Developed LabVIEW interface for the Rapid MPC prototyping: cRIO-9082 target deployement VI. is also used to easy tune the prediction and the control horizons and extract some needed target information. The last tab in this MPC controller VI indiates the establishing communication with the host PC and the successfully downloaded the cosimulation profile. It shows also the target IP address which be used inside the host VI to create a LabVIEW Network Stream data. After setting the IP address, the co-simulation profile is configured by coupling a setpoint to its applying time.

The second VI of the developed PIL platform, as depicted in Fig. 9, is built to display the tuned responses of the cRIO-9082 based controlled system in the host PC. Theses hardware Fig. 9: Developed LabVIEW interface for the Rapid MPC prototyping: plant hardware responses in the host PC. co-simulation results (plant output and control action) show the effectiveness of the proposed PIL methodology for rapid MPC prototyping. A 3D model of the MAGLEV system is also built and integrated within this VI in order to show, in the real-time, the suspended sphere motion.

On the other hands, setpoint, plant output and control action signals of the hardware MPC co-design are shown in Fig. 10 and Fig. 11, respectively. All performances of the implemented MPC algorithm for MAGLEV system, in terms of tracking and robustness, are guaranteed. The specified constraints on the control action signal are also respected. In order to highlight the proposed PIL co-simulation ability for choice and tuning of the MPC parameters, Fig. 12 and Fig. 13 show, respectively, the effect of varying the N c control and N p prediction horizons parameters on the step response of the controlled system. Values are indicated in the legend of each figure. So, from these figures we find out the usual results of the MPC strategy in terms of trajectory tracking and control law behavior. In Fig. 13, a very large prediction horizon has a reaction early and thus leads to a stable system, but less fast. However, a too short one causes a control law that reacts late, and with important values leading to the instability of the controlled system. Fig. 12 shows that increasing the control horizon parameter further improves the system response, but at some point they no longer earn anything. The values N c = 4, 5 or 8 give almost the same result. Thus, with a choice of parameters 0,5 0 0,02 0,06 0,08 0,1 0,12 0,14 0,16 0,18 0,2 0,22 0,24 0,26 0,28 0,3 0,32 0,34 0,36 0,38 0,4 0,42 0,44 0,46 0,48 0,5 0 0,02 0,04 0,06 0,08 0,1 0,12 0,14 0,16 0,18 0,2 0,22 0,24 0,26 0,28 0,3 0,32 0,34 0,36 0,38 0,4 0,42 0,44 0,46 0,48 N c = 5 and N p = 20, the PMC approach allowed to obtain a satisfactory tuning.

Setpoint

Clearly, this proposed LabVIEW/cRIO-based PIL platform also offers, in addition to the validation of the MPC hardware and software implementation materials, a solution for the selection and parameter settings of the MPC algorithms.

V. CONCLUSION

In this paper, a new CAD methodology for PIL cosimulation and rapid prototyping of model predictive controllers have been proposed and successfully applied to the position control of a MAGLEV system. The proposed SW/HW solution is based on the NI cRIO-9082 platform, associated to the LabVIEW/CDSim environment, and a host PC. The obtained results for such a LabVIEW/cRIO-based co-design approach are promising and open to improvement. Using such a CAD methodology for rapid MPC prototyping, users haves only to develop the host PC VI for new plant and introduce the new state space model as well as the MPC parameters in the target VI. Forthcoming works deal with Hardware-In-the-Loop (HIL) co-simulation of the MPC strategy and the final implementation of the control laws on the real MAGLEV benchmark available in our laboratory. In addition, the coupling with an optimization method to tune the weighting matrices to get the best possible performances in accordance with the specifications is also investigated.

Fig. 1 :

 1 Fig. 1: Prediction and control horizons in MPC framework.

Fig. 2 :

 2 Fig. 2: Model predictive control structure.

Fig. 3 :

 3 Fig. 3: LabVIEW CD Create MPC Controller VI.

Fig. 4 :

 4 Fig. 4: LabVIEW CD Implement MPC Controller VI.

Fig. 5 :

 5 Fig. 5: Block diagram of the cRIO-9082 based PIL cosimulation methodology.

Fig. 6 :

 6 Fig. 6: Used NI multi-core cRIO-9082 platform.

Fig. 7 :

 7 Fig. 7: Feedback MAGLEV 33-006 system.

Fig. 10 :Fig. 11 :

 1011 Fig. 10: Controlled position of the MAGLEV suspended sphere.

Fig. 12 :

 12 Fig. 12: Effect of varying the control horizon N c for a fixed prediction horizon (N p = 20).

Fig. 13 :

 13 Fig. 13: Effect of varying the prediction horizon N p for a fixed control horizon (N c = 5).