Xabier Jaureguiberry

Emmanuel Vincent

Gaël Richard

Multiple-order Non-negative Matrix Factorization for Speech Enhancement

Keywords: Variational Bayes, Non-negative Matrix Factorization, Model Averaging, Speech Enhancement

Amongst the speech enhancement techniques, statistical models based on Non-negative Matrix Factorization (NMF) have received great attention. In a single channel configuration, NMF is used to describe the spectral content of both the speech and noise sources. As the number of components can have a crucial influence on separation quality, we here propose to investigate model order selection based on the variational Bayesian approximation to the marginal likelihood of models of different orders. To go further, we propose to use model averaging to combine several single-order NMFs and we show that a straightforward application of model averaging principles is inefficient as it turned out to be equivalent to model selection. We thus introduce a parameter to control the entropy of the model order distribution which makes the averaging effective. We also show that our probabilistic model nicely extends to a multipleorder NMF model where several NMFs are jointly estimated and averaged. Experiments are conducted on real data from the CHiME challenge and give an interesting insight on the entropic parameter and model order priors. Separation results are also promising as model averaging outperforms single-order model selection. Finally, our multiple-order NMF shows an interesting gain in computation time.

Introduction

Speech enhancement has received great attention since it is at stake in numerous industrial applications. The literature proposes a variety of methods which fall into two categories [START_REF] Benesty | Speech enhancement[END_REF] : multichannel vs. single-channel. In a single-channel configuration, which is the focus of this paper, spatial diversity cannot be exploited. Single-channel techniques can be classified into three main categories [START_REF] Loizou | Speech enhancement: theory and practice[END_REF]: spectral subtractive algorithms [START_REF] Boll | Suppression of acoustic noise in speech using spectral subtraction[END_REF], subspace algorithms [START_REF] Ephraim | A signal subspace approach for speech enhancement[END_REF] and statistical-model-based algorithms [START_REF] Ephraim | Statistical-model-based speech enhancement systems[END_REF]. Amongst the latter, various models aim at describing the spectral content of both the speech and background noise sources in order to be able to distinguish between them. To build such models, one can resort for instance to Gaussian mixture models [START_REF] Hao | Speech enhancement, gain, and noise spectrum adaptation using approximate Bayesian estimation[END_REF], exemplar-based methods [START_REF] Gemmeke | HMMregularization for NMF-based noise robust ASR[END_REF] or codebook-driven techniques [START_REF] Srinivasan | Codebook driven short-term predictor parameter estimation for speech enhancement[END_REF][START_REF] Mowlaee | A joint approach for single-channel speaker identification and speech separation[END_REF]. Non-negative Matrix Factorization (NMF) is one of the most popular class of source models [START_REF] Raj | Phoneme-dependent NMF for speech enhancement in monaural mixtures[END_REF][START_REF] Mohammadiha | Single channel speech enhancement using Bayesian NMF with recursive temporal updates of prior distributions[END_REF][START_REF] Moritz | Noise robust distant automatic speech recognition utilizing NMF based source separation and auditory feature extraction[END_REF][START_REF] Geiger | The TUM + TUT + KUL approach to the 2nd CHiME challenge: Multi-stream ASR exploiting BLSTM networks and sparse NMF[END_REF] and it has achieved great performance in the latest CHiME contest [START_REF] Geiger | The TUM + TUT + KUL approach to the 2nd CHiME challenge: Multi-stream ASR exploiting BLSTM networks and sparse NMF[END_REF][START_REF] Vincent | The second CHiME speech separation and recognition challenge: An overview of challenge systems and outcomes[END_REF].

The order of an NMF model, also called the number of components, is known to have a noticeable influence on separation quality [START_REF] Bertin | Blind signal decompositions for automatic transcription of polyphonic music: NMF and K-SVD on the benchmark[END_REF]. However, there is a few literature about how to determine the best number of components [START_REF] Mowlaee | A joint approach for single-channel speaker identification and speech separation[END_REF] and choosing it is often driven by an experimental assessment. The introduction of a statistical formulation of NMF [START_REF] Virtanen | Bayesian extensions to non-negative matrix factorisation for audio signal modelling[END_REF] has made the applica-tion of model selection principles possible. To do so, the literature advocates the use of a full Bayesian framework [START_REF] Attias | A variational Bayesian framework for graphical models[END_REF][START_REF] Hoeting | Bayesian model averaging: a tutorial[END_REF] in order to compute the marginal likelihood of a model, also named the evidence. For a given task, once the marginal likelihood has been estimated for several models of different orders, the model with the highest marginal likelihood is reputed to be the most likely model to explain the observation. In practice, a full Bayesian treatment of NMF is intractable and approximate inference is required instead. In particular, Variational Bayesian (VB) inference is becoming popular since it is less computationally demanding than sampling methods [START_REF] Attias | A variational Bayesian framework for graphical models[END_REF][START_REF] Cemgil | Bayesian inference for nonnegative matrix factorisation models[END_REF].

VB proposes to approximate the marginal likelihood by a lower bound called the free energy. This free energy can then be used in place of the true marginal likelihood for model order selection. Two types of approaches have been exploited so far: parametric methods which consist in computing several NMFs and selecting the one which has the highest free energy [START_REF] Cemgil | Bayesian inference for nonnegative matrix factorisation models[END_REF], and nonparametric methods which consider a single NMF with a potentially infinite number of components and which iteratively deactivates the irrelevant components [START_REF] Hoffman | Bayesian nonparametric matrix factorization for recorded music[END_REF].

To go further, we propose here to apply model averaging principles to NMF models [START_REF] Hoeting | Bayesian model averaging: a tutorial[END_REF]. Indeed, the study in [START_REF] Jaureguiberry | Introducing a simple fusion framework for audio source separation[END_REF] has shown that it is worth combining several NMFs of different orders instead of selecting a unique one. The free energy given by VB inference can be used to compute the posterior probability of each number of components and to weight each NMF [START_REF] Attias | A variational Bayesian framework for graphical models[END_REF]. We also derive from our averaging of single-order NMFs a novel multiple-order model which jointly estimates and averages several NMFs of different orders. However, our contribution underlines that a straightforward application of model averaging based on the free energy is inefficient as it turns out to select a single NMF. To avoid this behaviour, we propose to use a parameter which controls the entropy of the distribution of the number of components. As this entropic parameter and the order priors need to be chosen beforehand, we propose to learn them on a training database. Our experiments conducted on real data from the CHiME challenge [START_REF] Vincent | The second 'CHiME' speech separation and recognition challenge: datasets, tasks and baselines[END_REF] show promising results as both the multiple-order NMF and the averaging of singleorder NMFs outperform single-order model selection, thanks to the introduction of the entropic parameter. Moreover, our multiple-order NMF turns out to be less computationally expensive than the averaging of single-order models for equivalent performance.

In the following, Section 2 will be dedicated to the presentation of the single-order NMF model and its averaging. In Section 3, we will introduce our novel multiple-order NMF model. The entropic parameter will be presented in Section 4 before being experimentally evaluated in Section 5. Section 6 will give a concise conclusion.

Single-order NMF

Our single-order NMF model is a single-channel simplified formulation of the model exposed in [START_REF] Adilȏglu | Variational Bayesian inference for source separation and robust feature extraction[END_REF]. The degraded speech signal is supposed to be a linear mixture of a clean speech sig-nal and a background noise signal. As such, the mixing equation can be written in the short-time Fourier transform (STFT) domain as

x f n = As f n + f n = s 1,f n + s 2,f n + f n (1)
in which f denotes the frequency bin and n the time frame, s 1,f n is the target speech signal, s 2,f n is the background noise and f n represents sensor noise. By denoting the source vector as s f n = [s 1,f n s 2,f n ] T , the mixing equation can also be written in a matrix form thanks to the mixing matrix A = [1 1]. Such a formulation permits future extension of our model to more sources and channels.

Probabilistic model

Speech and background noise sources are both treated in the same way and indexed by j = {1, 2}. We assume that the source s j,f n follows a circularly-symmetric complex normal distribution s j,f n ∼ N (0, v j,f n ) whose variance is the result of an NMF so that

v j,f n = K j k=1 w j,f k h j,kn . (2)
Kj is the order of the NMF, a.k.a. the number of components.

In this single-order NMF formulation, it is an hyperparameter to be chosen. The variance of source j can also be written in a matrix form so that Vj = WjHj. Wj and Hj are commonly called the dictionary and the activation matrix. To achieve full Bayesian inference, the parameters of both the dictionary and activation matrix are seen as random variables as well. Contrary to [START_REF] Adilȏglu | Variational Bayesian inference for source separation and robust feature extraction[END_REF], we assume as in [START_REF] Hoffman | Bayesian nonparametric matrix factorization for recorded music[END_REF] that the NMF parameters follow a Gamma distribution

w j,f k ∼ Γ(a, a) , h j,kn ∼ Γ(b, b) (3)
where a and b are hyperparameters to be chosen. Finally, sensor noise is supposed to follow a Gaussian distribution of variance σ 2 . Denoting X = {x f n } n=1..N f =1..F and S = {s f n } n=1..N f =1..F for the sake of readability, the log-likelihood can be formulated as

log p(X|S) = N n=1 F f =1 log N x f n |As f n , σ 2 .
(4)

Variational inference

In the following, we denote the set of all model parameters as Z = {S, W, H} with W = {Wj} j=1,2 and H = {Hj} j=1,2 . VB inference aims at approximating the posterior distribution of the model parameters p(Z|X) with a factorized variational distribution q(Z) which is here defined as q(Z) = q(S)q(W)q(H) so that

q(Z) = f n q(s f n ) j,f k q(w j,f k ) j,kn q(h j,kn ). (5)
In VB inference, minimizing the Kullback-Leibler divergence between q(Z) and p(Z|X) is equivalent to approximating the log marginal likelihood log p(X) by the so-called free energy L [q] defined as

L [q] = q(Z) log p(X, Z) q(Z) dZ = E Z log p(X, Z) q(Z) (6)
with p(X, Z) = p(X|S) p(S|W, H) p(W) p(H) being the joint distribution. VB inference then consists in iteratively maximizing the free energy with respect to each factor in [START_REF] Ephraim | Statistical-model-based speech enhancement systems[END_REF]. In practice, the computation of the free energy is intractable and it needs to be further approximated by a parametric lower bound B [q]. For a detailed explanation, the reader is referred to [START_REF] Adilȏglu | Variational Bayesian inference for source separation and robust feature extraction[END_REF].

Deriving B [q] with respect to each factor leads to the update rules. The variational distribution of the sources is identified to a bivariate Gaussian distribution [START_REF] Adilȏglu | Variational Bayesian inference for source separation and robust feature extraction[END_REF]

q(s f n ) = N µ s,f n , Σ s,f n with parameters µ s,f n = Σ s,f n A T x f n σ 2 , Σ s,f n = C -1 f n + 1 σ 2 J -1 (7)
where J is a matrix of ones of size 2 × 2,

C f n = diag(C j,f n )j=1,2 and C j,f n = K j k=1 E Z\s j 1 w j,f k h j,kn -1 (8)
The notation E Z\Z i [.] denotes the expectation over all model parameters Z except Zi.

The variational distributions of the NMF parameters are identified to generalized inverse Gaussian (GIG) distributions [START_REF] Hoffman | Bayesian nonparametric matrix factorization for recorded music[END_REF] which are controlled by three parameters τ , ρ and γ. The updates of these parameters for the activation matrix Hj of source j are given by:

τ H j = E 1 Hj .2 •   E 1 Wj .-1 T E |Sj| .2 • C .-2 j   ρ H j = b + E [Wj] T E [Vj] .-1 , γ H j = b (9)
where the notation • denotes the Hadamard product, M .x and M T respectively denote element-wise exponentiation and transposition of matrix M. Cj is the matrix composed of the coefficients C j,f n defined in [START_REF] Srinivasan | Codebook driven short-term predictor parameter estimation for speech enhancement[END_REF]. Note that the same update rules can be found for Wj by replacing and reordering the terms accordingly.

In a speech enhancement context, the variable in which we are interested is the clean speech source s 1,f n . The STFT coefficients of the estimated sources are given by the mean µ s,f n of the posterior distribution q(s f n ) in [START_REF] Gemmeke | HMMregularization for NMF-based noise robust ASR[END_REF]. For the speech source s 1,f n , this expectation simplifies to

µs 1,f n = C 1,f n C 1,f n + C 2,f n + σ 2 x f n . ( 10
)
We recognize the classical expression of the Wiener filter where deterministic estimates of the source power spectra have been replaced by the expectations C j,f n .

Model averaging

We now assume that the above single-order NMF framework has been used with M different models. The model m is defined by its order denoted Km = {K1m, K2m}, in which K1m (resp. K2m) is the number of components of the speech source (resp. background noise source). The posterior distribution qm(Z) has thus been estimated for each model m = 1..M . Bayesian model averaging [START_REF] Hoeting | Bayesian model averaging: a tutorial[END_REF] proposes to average these posterior probabilities with respect to the posterior probability p(Km|X) of each model. Thanks to Bayes' rule, this posterior probability can be expressed as the product of the prior probability πm of model m and its likelihood p(X|Km) so that p(Km|X) ∝ πm p(X|Km).

(

) 11
As we have already highlighted, the computation of the likelihood p(X|Km) is intractable. However, the choice of a VB framework gives us the opportunity to replace it with the free energy expressed in (6) so that p(X|Km) ≈ exp(Lm).

Applying model averaging to the M posterior distributions of the sources leads to the new estimate

q(s f n ) = 1 δ M m=1 πme Lm qm(s f n ) (12)
where qm(s f n ) is the source posterior distribution estimated for model m and δ = M m=1 πme Lm aims at normalizing the posterior probability of Km so that it sums to 1. The STFT coefficients of the estimated speech source is now given by a linear combination of the expectations µs 1m,f n computed for each model m as in [START_REF] Raj | Phoneme-dependent NMF for speech enhancement in monaural mixtures[END_REF]. This formulation is equivalent to the temporal fusion by linear combination introduced in [START_REF] Jaureguiberry | Introducing a simple fusion framework for audio source separation[END_REF] but differs from it in that the fusion coefficients now depend on the signal to be processed through the free energy Lm.

Multiple-order NMF

Our novel multiple-order NMF model introduced here aims at jointly estimating and averaging several NMFs of different orders. For more details on the model and the derivation of the VB inference, the reader is referred to [START_REF] Jaureguiberry | Variational Bayesian model averaging for audio source separation[END_REF].

Probabilistic model

Contrary to the single-order NMF model of Section 2, the number of components Kj of source j is now seen as a random variable which follows a categorical distribution Kj ∼ Cat πj1, ..., πjm, ..., πjM j [START_REF] Geiger | The TUM + TUT + KUL approach to the 2nd CHiME challenge: Multi-stream ASR exploiting BLSTM networks and sparse NMF[END_REF] where m indexes the Mj possible number of components Kj1, ..., Kjm, ..., KjM j , each having an a priori probability of πjm. As above, we assume that each number of components has its specific NMF parameters so that in the following, they will be indexed with m such as w jm,f k and h jm,kn . The priors on the NMF parameters and the sources remain unchanged in comparison with Section 2. The main difference is that a single posterior distribution of the sources is now estimated.

By denoting

K = {K1, K2}, W = {Wjm} m=1..M j j=1,2 , H = {Hjm} m=1..M j j=1,2
and p(K) = p(K1)p(K2), the joint distribution becomes p(X, Z, K) = p(X|S) p(S|W, H) p(W|K) p(H|K) p(K).

Variational inference

In order to include the variables related to the numbers of components, the variational distribution of ( 5) is modified as follows q(Z, K) = q(S) q(W|K) q(H|K) q(K).

(14) By minimizing the corresponding free energy, the posterior probability of the number of components q(Kj) is obtained as

log q(Kj) = E Z\K j [log p(X, Z, K)] -E Z\K j [log q(Wj|Kj)] -E Z\K j [log q(Hj|Kj)] + const.
By developping, reordering and taking the exponential of these terms, the posterior probability of Kj can be formulated as ∀m, q(Kjm) ∝ πjm e L jm [START_REF] Bertin | Blind signal decompositions for automatic transcription of polyphonic music: NMF and K-SVD on the benchmark[END_REF] in which the term Ljm = E Z\K j log p(X,Z) q(Z)

is similar to the free energy expressed in [START_REF] Hao | Speech enhancement, gain, and noise spectrum adaptation using approximate Bayesian estimation[END_REF] in the single-order NMF case. The inference of the variational distributions over the NMF parameters remains unchanged in comparison to the single-order NMF model and the variational distribution of the sources is still identified as a bivariate Gaussian distribution with the parameters defined in [START_REF] Gemmeke | HMMregularization for NMF-based noise robust ASR[END_REF]. The only difference relies in the term C -1 f n which is now a linear combination of terms related to each NMF order so that, with δ =

M j m=1 πjme L jm , ∀j, C -1 j,f n = 1 δ M j m=1 πjme L jm C -1 jm,f n . ( 16
)

Controlling the order posterior entropy

We introduced two NMF frameworks which both combine several NMFs of different orders. In Section 2, the scheme consists in computing several single-order NMFs and combining them after the variational inference thanks to the linear combination in [START_REF] Moritz | Noise robust distant automatic speech recognition utilizing NMF based source separation and auditory feature extraction[END_REF]. In Section 3, the scheme jointly estimates and average several NMFs according to [START_REF] Virtanen | Bayesian extensions to non-negative matrix factorisation for audio signal modelling[END_REF]. However, the averaging is similar in both cases as the averaging weights of ( 12) and ( 16) are of the form πm exp(Lm).

To enforce this comparison, we now assume as a particular case that the number of components of the background noise K2 is fixed. In the single-order NMF framework, the averaging rule [START_REF] Moritz | Noise robust distant automatic speech recognition utilizing NMF based source separation and auditory feature extraction[END_REF] remains unchanged. The model m is now entirely defined by the number of components K1m of the speech source. In the multiple-order case, the number of components K1 of the speech source is now the only one to be considered as a random variable and the average is only effective on the related NMF parameters and not on the background noise side. Hence, ( 16) only holds for j = 1 and by dropping the index j when unnecessary, it becomes

C -1 1,f n = 1 δ M m=1 πme Lm C -1 1m,f n .
As a consequence, Lm and πm now have the exact same signification in both cases. πm denotes the prior probability of model m of order K1m whereas Lm denotes its free energy. Note however that in the multiple-order case, all speech models share the same background noise estimate whereas in the single-order case, one background noise model is estimated per order K1m.

As stated in [START_REF] Jaureguiberry | Introducing a simple fusion framework for audio source separation[END_REF][START_REF] Jaureguiberry | Variational Bayesian model averaging for audio source separation[END_REF], it is worth combining NMFs of different orders instead of selecting a unique order as it can improve the separation performance. However, our models, which apply model averaging in a straightforward way, are not achieving that goal. Indeed, preliminary tests have shown that the free energies Lm being very large valued, the posterior probabilities q(K1m) are all equal to zero except the one which depends on the highest free energy and which is thus equal to 1. As such, model averaging turns out to be model selection. To avoid this behaviour, we propose to scale the free energies Lm by a factor β ≥ 1 to be determined. Hence, the posterior probability of order K1m becomes q(K1m) ∝ πm e Lm/β . (17) This is equivalent to controlling the entropy of the distribution q(K1m) in a way similar to [START_REF] Katahira | Deterministic annealing variant of variational Bayes method[END_REF]. Small values of β will favor peaky distributions with one q(K1m) close to 1 as in model selection, whereas higher values of β will result in a more uniform distribution. The determination of the values of β and πm will be addressed in Section 5.

Experiments

In this section, we propose to evaluate and compare both the multiple-order NMF model and the model averaging of singleorder NMFs. We also explore oracle and learning methods to determine the prior probabilities of the number of components πm as well as the entropic parameter β. To do so, we rely on the PASCAL CHiME corpus [START_REF] Vincent | The second 'CHiME' speech separation and recognition challenge: datasets, tasks and baselines[END_REF] which features recordings of real domestic noise and speech utterances from diverse speakers.

Dataset

The CHiME corpus is divided in three datasets using 34 distinct speakers. Firstly, the training set is composed of 500 utterances of each speaker in reverberated conditions. This allows us to learn an NMF dictionary W1m to describe each speaker by applying a single-order NMF on the concatenation of these utterances with a chosen number of components K1m. All these learned dictionaries hence describe the same source but at different levels of details. The other datasets, namely the develop- Table 1:Average SDR (dB) and computation time for VB and oracle selections, simple-order NMF averaging and multiple-order NMF ment and test sets, are both composed of 600 reverberated utterances. Each utterance, pronounced by one of the 34 speakers, has been mixed within a noise background so that the mixture of clean speech and noise achieved a given signal-to-noise ratio (SNR) amongst: -6, -3, 0, 3, 6 and 9 decibels (dB). The noise background signals come from real recordings of a domestic living room and the utterances have been placed at controlled time stamps according to the desired SNR. This allows us for each utterance to select 10 seconds of noise without speech, after and/or before the utterance, in order to learn an NMF dictionary W2 in a way similar to the learning of the NMF speaker models. We chose a single-order with a fixed number of components K2 = 16.

In the following, we have used the development set for learning purpose, whereas for evaluation, we have randomly selected in the test set a total of 24 utterances, i.e., 4 by SNR. As the original data were two channel mixtures, we will work here on the mean of both channels to restrain our study to the single channel case. Both the single-order NMFs and the multiple-order NMF have been used with numbers of components K1m = 2 m with m = 1..7 for the speaker source. The dictionaries of noise and speech NMFs have been fixed to their learned values and the activation matrices have been initialized with the mean activation values estimated in the corresponding learning step. The hyperparameter of the Gamma prior related to the activation matrices has been fixed to b = 0.2 and sensor noise is assumed to be of variance σ 2 = 10 -6 . Finally, the separation quality has been evaluated by the signal-to-distortion ratio (SDR) of the target speech, expressed in decibels [START_REF] Vincent | Performance measurement in blind audio source separation[END_REF].

Learning the priors and the entropic parameter

In order to average the NMFs, we need to determine the prior probabilities πm of each number of components as well as the entropic parameter β. To do so, we propose two approaches.

In particular, we propose to learn both πm and β thanks to the development set. Indeed, as the original reverberated signals are available, we can evaluate the performance of the single-order NMF scheme on each mixture l of the development dataset and for each number of components K1m = 2 m . We can thus find the set of πm and β which maximizes the mean SDR of the speech signals over the L examples of the development set. This is equivalent to solving the non-linear optimization problem

argmin (π 1 ,...,π 7 ,β) l,f n µ (l) s 1,f n - s(l) 1,f n 2 (18)
in which s(l) 1,f n is the original reverberated speech signal of example l, µ

(l) s 1,f n = m πm exp(Lm/β)µ (l) s 1m,f n and µ (l) s 1m,f n
is defined for example l as in [START_REF] Raj | Phoneme-dependent NMF for speech enhancement in monaural mixtures[END_REF] for the number of components K1m. Note that such a learning is impossible in the multipleorder NMF case as the averaging rule ( 16) is computed at each iteration of the inference algorithm. However, these learned values of πm and β can be used in both the averaging of singleorder NMFs and the multiple-order NMF model. For computational convenience, we have restrained the size of the development set for learning. We thus propose to learn the parameters on L = 36 randomly picked examples, i.e., 6 by SNR.

In order to evaluate the performance of the learned πm and β, we can also solve the optimization problem [START_REF] Hoeting | Bayesian model averaging: a tutorial[END_REF] for the example of the test set being processed. However, it is worth noting that such a result is unreachable in practice and is thus denoted as the oracle result. As in the learning case, oracle πm and β will be used in both the single-order and multiple-order schemes.

Results and comments

Table 1 shows the results of our study grouped in three categories. The SDRs are averaged for each SNR over the 4 selected examples as well as over the whole selected test set. The baseline is given by the VB selection result, which consists in selecting for each example the single-order model that has the highest free energy Lm, and the oracle selection result which, on the basis of the original reverberated speech signal, consists in selecting for each example the single-order model which gives the best SDR. Note that VB selection is equivalent to the averaging of single-order NMFs without using the entropic parameter we proposed in Section 4. The second group of results is based on the oracle πm and β values whereas the third group is based on the learned πm and β values, as explained in Section 5.2. For each example of the test set, these oracle and learned parameters have been used for both the single-order NMF (so-NMF) averaging and the multiple-order NMF (mo-NMF). These results first show that regardless of the SNR, the VB selection always fails to select the best model in term of SDR. In average, the VB selection underperforms oracle selection by 1.3 dB. Oracle results show that the averaging of single-order NMFs can always outperform both VB and oracle selection. The multiple-order NMF using these same oracle parameters implies a loss of only 0.35 dB in comparison with the singleorder averaging. However, this loss is somehow counterbalanced by an interesting gain in computation time as mo-NMF is 4 times faster than the averaging of so-NMFs.

Recalling that oracle methods are not reachable in practice, the results based on learned πm and β show an interesting performance. Indeed, both single-order and multiple-order methods outperform VB selection which is the only other method to be practicable. Moreover, both methods nearly reach the SDR of the oracle selection. Our proposed method to learn πm and β is thus efficient as it allows to improve separation performance by almost 1.1 dB. We can also notice that in realistic conditions, the multiple-order NMF gives similar SDRs to the averaging of single-order NMFs but it is far less computationally demanding.

Conclusion

We introduced two different NMF frameworks which aim at combining several NMFs of different orders: a single-order model scheme which applies Bayesian model averaging and a novel multiple-order model which jointly estimates and averages several NMFs. The parameter we introduced to control the entropy of the order posterior has been shown to be essential in order to make the averaging effective in both models. Finally, experimental results on real data showed that both models outperform traditional VB selection, the multiple-order NMF being furthermore less computationally demanding. Future works will focus on the study of a time-varying model averaging.

This work was partly supported under the research programme EDiSon3D (ANR-13-CORD-0008-01) funded by ANR, the French State agency for research.


Olivier Dorival
email: dorival@lmt.ens-cachan.fr

Philippe Rouch
email: rouch@lmt.ens-cachan.fr

Olivier Allix
email: allix@lmt.ens-cachan.fr

A SUBSTRUCTURED VERSION OF THE VARIATIONAL THEORY OF COMPLEX RAYS DEDICATED TO THE CALCULATION OF ASSEMBLIES WITH DISSIPATIVE JOINTS IN THE MEDIUM-FREQUENCY RANGE

Keywords: Joint Modeling, Medium Frequency, Damping, VTCR, Substructuration, Variational Approach

A substructured version of the variational theory of complex rays dedicated to the calculation of assemblies with dissipative joints in the medium-frequency range

Introduction

The prediction of the response of a complex structure in the medium-frequency (MF) range presents two major difficulties. The first difficulty is that the solution of a problem involves small-wavelength phenomena. Therefore, using the finite element method (FEM), the calculation of this solution leads to a huge number of degrees of freedom. Indeed, at least seven elements per wavelength are necessary to represent oscillating solutions properly. Moreover, since k 3 h 2 must remain constant [START_REF] Deraemaeker | Dispersion and pollution of the FEM solution for the Helmholtz equation in one, two and three dimensions[END_REF], the cost of FE analysis increases with the frequency. The second difficulty is that due to the high modal density in the MF domain the response is governed mainly by the level and the nature of the dissipative phenomena, whose mechanisms are very difficult to describe, model and measure. In the case of thin structures -such as car chassis, satellites or submarines -, damping is localized mainly in the joints between substructures, which can be of different types.

Many research works have been devoted to answering the first point. Some aim at extending the classical low-frequency approaches (FEM) to the medium-frequency range [START_REF] Farhat | The Discontinuous Enrichment Method[END_REF][START_REF] Farhat | A Discontinuous Galerkin Method with Lagrange Multipliers for the Solution of Helmholtz Problems in the Mid-Frequency Regime[END_REF][START_REF] Harari | Improved Finite Element Methods for Elastic Waves[END_REF]. Other approaches, such as Statistical Energy Analysis (SEA) [START_REF] Lyon | Power flow between linearly coupled oscillators[END_REF], try to relax some of the assumptions associated with high-frequency methods in order to become applicable to the MF domain.

In this study, we use the Variational Theory of Complex Rays (VTCR) [START_REF] Ladevèze | The variational theory of complex rays for the calculation of mediumfrequency vibrations[END_REF], which deals with the first difficulty by using two-scale shape functions which verify the dynamic equations and the constitutive law of the substructures. This leads to a small, frequency-independent linear system [START_REF] Rouch | The variational theory of complex rays : a predictive tool for medium-frequency vibrations[END_REF][START_REF] Ladevèze | Calculation of Medium-Frequency Vibrations over a Wide Frequency Range[END_REF]. We focus particularly on the second difficulty.

Noting that one of the specificities of the VTCR is its variational formulation of all the boundary conditions between substructures, we concentrate on the modeling of joints between substructures: this feature allows us to develop the equations of complex interfaces involving a dissipative part easily, which is necessary to obtain accurate predictions of the behavior of real structures.

To improve this process, we derive a substructured version of the VTCR which allows us to treat interfaces with greater flexibility. Indeed, in the classical version of the VTCR, all the interface conditions of the entire structure are coupled, which can be exceedingly involved to improve only some joint equations, for example when updating the transmission conditions of the structure using experimental data or when analyzing multiresolution problems. This substructured reformulation is based on a condensation approach which consists in separating the properly modeled substructures and boundaries from the more inaccurate complex joints. After the introduction of interface degrees of freedom, one ends up with a constrained minimization problem.

Finally, in order to predict the behavior of real structures accurately, we introduce models of complex joints involving a dissipative part. These models were implemented into the VTCR. We present applications concerning 3D plate assemblies connected by heterogeneous joints with distributed parameters (such as mass, stiffness and damping) which can be identified on generic assemblies.

A substructured version of the VTCR

Introduction

The classical version of the VTCR [START_REF] Ladevèze | The variational theory of complex rays for the calculation of mediumfrequency vibrations[END_REF]), which will not be detailed in this paper, is entirely dedicated to the modeling of joints, since it is based on the decomposition of the structure into simple substructures, which involves transmission conditions. However, this formulation treats all the boundaries and interfaces of the structure globally, which makes it cumbersome for studying the vibrations of the same structure with different joint equations if, for example, one wants to model the interfaces using stochastic interface parameters to describe the dispersion due to the production process. Therefore, after introducing the reference problem, we derive a substructured version of the VTCR.

Many substructuration methods have appeared from the beginning of numerical methods. The first studies were based on modal analysis of each substructures (see the Com-ponent Modal Synthesis (CMS) [START_REF] Hurty | Dynamic Analysis of Structural Systems Using Component Modes[END_REF]). The Craig-Bampton approach proposes to compute substructures modes with clamped interfaces [START_REF] Craig | Coupling of Substructures for Dynamic Analysis[END_REF]). In the Mc-Neal method, the substructures modes are computed with free interface conditions [START_REF] Mcneal | A Hybrid Method For Component Mode Synthesis[END_REF]. Farhat and Géradin proposed a method allowing non compatible interfaces [START_REF] Farhat | On a Component Mode Synthesis Method and its Application to Incompatible Substructures[END_REF].

Other domain decomposition methods are not based on modal analysis, but introduce various interface quantities : the Balancing Domain Decomposition method [START_REF] Roux | Méthodes de Résolution Par Sous-Domaines En Statique[END_REF][START_REF] Mandel | Balancing Domain Decomposition[END_REF]) uses primal quantities, whereas the FETI method [START_REF] Farhat | A Method for Finite Element Tearing and Interconnecting and its Parallel Solution Algorithm[END_REF], Farhat, et al. 1994[START_REF] Farhat | The Two-Level FETI Method for Static and Dynamic Plate Problems Part I: An Optimal Iterative Solver for Biharmonic Systems[END_REF]) uses dual interface quantities. Hybrid methods like the LATIN approach [START_REF] Ladevèze | Nonlinear Computational Structural Method : New Approaches and Non-Incremental Methods of Calculations[END_REF] propose not to favour dual or primal quantities.

The reference problem

Here, for the sake of simplicity, the problem will be formulated for an assembly of only two substructures, but this can be easily generalized to an assembly of n substructures. As an illustration, let us consider two thin, homogeneous, isotropic and elastic Kirchhoff-Love plates and study the steady-state vibrations of this assembly. Classically, all quantities are defined in the complex domain and depend on the space quantities alone: the complex quantity Q (X) is associated with Q (X, t) by the relation Q (X, t) = Q (X) . exp(i ω t), where ω is the fixed angular frequency.

S 1 w 1, M 1  12 S 2 w 2, M 2 K d M d
Fig. 1. The reference problem

Let S 1 and S 2 be two plates, ∂S 1 and ∂S 2 their respective boundaries solicited harmonically at a fixed angular frequency ω, and Γ 12 the interface between S 1 and S 2 , which, for clarity, will first be considered to be perfectly rigid. The boundary conditions are the following, with l = 1, 2: prescribed out-of-plane displacement w d l on a part ∂ w d S l of ∂S l , prescribed slope w d n, l on ∂ w d n S l , prescribed bending moment M d l on ∂ M d S l and, finally, prescribed Kirchhoff shear K d l on ∂ K d S l (see Fig. 1). The quantities of interest are the displacements w l and the moments M l .

For S l , let us introduce the admissibility space S l ad of the displacement-moment pairs s l = (w l , M l ) defined on S l such that:

s l ∈ S l ad ⇔                        w l ∈ U (the set of the finite-energy displacement fields H 2 (S) )
M l ∈ S (the set of the finite-energy moment fields [L 2 (S) ] 3 )

∆∆w l -k 4 l w l = 0 on S l M l = 2h 3 l 3 (1 + iη l ) K l, P S X (w l ) for l = 1, 2 (1)
with

k 4 l = 3ρ l ω 2 (1 -ν 2 l ) (1 + iη l )E l h 2 l
where k l is the wave number and ρ l , ν l , η l , E l and 2h l denote respectively the density, Poisson's ratio, structural damping coefficient, Young's modulus and thickness of Plate S l . K l, P S is Hooke's tensor for plane stress and X is the curvature operator. The fields of S l ad are admissible in the sense that they verify both the local equilibrium and the constitutive relation exactly. The reference problem is:

Problem 1. Find s l such that: s l ∈ S l ad              w l = w d l on ∂ w d S l w l, n l = w d n l on ∂ w d n S l M l = M d l on ∂ M d S l K l = K d l on ∂ K d S l w 1 = w 2 on Γ 12 w 1, n 1 = w 2, n 2 on Γ 12 M 1 + M 2 = 0 on Γ 12 K 1 + K 2 = 0 on Γ 12 for l=1, 2 (2) with M l = n l M l n l and K l = n l div[M l ] + (t l M l n l ) ,t l
The reference problem has one solution if [START_REF] Ladevèze | The variational theory of complex rays for the calculation of mediumfrequency vibrations[END_REF]):

• Hooke's tensors are positive definite, • η l > 0.

A substructured version of the VTCR

Let us introduce a substructured version of the VTCR which consists in separating the properly modeled substructures from the less accurate complex joints. The cornerstone of this new formulation is the introduction of interface quantities, which we denote •. Then, the interface equations can be solved by minimizing a functional under a set of constraints expressing the compatibility between substructures and interfaces in the form of classical VTCR problems, which will be briefly reminded. In doing this, we ensure that the solution within the substructures is exact: the only error, which is perfectly acceptable, is in the VTCR discretization of the substructures (anterior studies [START_REF] Rouch | The variational theory of complex rays : a predictive tool for medium-frequency vibrations[END_REF] have focused on the influence of the VTCR discretization and it will not be the subject of the present paper). Besides, the joints are dealt with less accurately through the minimization problem. This reformulation gives us greater flexibility to modify the interface equations: the substructures' minimization problem with constraints is built once and for all, and the calculation with a different set of joints requires changes only in the minimized interface residual. In future works, we will take advantage of this new formulation to calculate stochastic joint behavior, to deal with multiresolution problems or to identify joint parameters. In this paper, we concentrate on the presentation of this new formulation.

Formulation of the continuous substructured VTCR

This reformulation is based on the separation of the equations of the joints through the introduction of interface quantities. For clarity, we present only the separation of the interface's bending behavior. Considering the interface Γ 12 , we introduce the interface quantities (w 1, n 1 , M 1 ) and (w 2, n 2 , M 2 ) of Substructures 1 and 2 respectively. Each interface quantity (w l, n l , M l ) must be equal to the corresponding substructure quantity (w l, n l , M l ) (at least in the continuous case). The interface, which can be viewed as a real entity, has its own behavior described by an equilibrium equation and a constitutive law involving the interface quantities (w 1, n 1 , M 1 ) and (w 2, n 2 , M 2 ). Thus, the basic problem is divided into two subproblems (see Fig. 2): Substructure l must verify a substructure/interface compatibility problem while the interface quantities must minimize a L 2 residual of joint equations.

The substructure/interface compatibility problem

Among several possible options, we chose a condensation approach with which prescribed interface quantities M l :

 12 w 2 w 1 M 1 M 2 S 2 w 2, M 2 K d M d M 2 w 2 S 1 w 1, M 1 M 1 w 1

Fig. 2. The VTCR problem with substructuration

Substructure l/interface compatibility problem:

-dynamic equation -constitutive relation -boundary conditions: prescribed slope:

w d n, l = w d n, l on ∂ w d n S l prescribed bending moments: M l = M d l on ∂ d M S l prescribed interface quantities: M l = M l on Γ 12 w n, l = w n, l on Γ 12 (3)
This problem is solved through the use of the VTCR, which we are about to introduce. The VTCR is primarily a global formulation of the boundary conditions and transmission conditions in terms of both displacements and moments:

Problem 2 l . Find s l in S l
ad such that:

A ∂S l (δs l , s l ) -L ∂S (δs l ) + C Γ l δs l , s l , s l = 0 ∀δs l ∈ S l ad for l = 1, 2, (4)
where:

A ∂S l (δs l , s l ) -L ∂S l (δs l ) = R e -iω - ∂ w d S l δK l (w l -w d l ) * dL + ∂ w d n S l δM l (w l, n l -w d l, n l ) * dL + ∂ M d S l (M l -M d l ) δw * l, n l dL - ∂ K d S l (K l -K d l ) δw * l dL (5)
and with the prescribed interface quantities on Γ:

C Γ = R e -iω Γ 12 (M l -M l ) δw * l, n l dL (6)
R e (Q) and Q * designate the real part of Q and the conjugate of Q respectively. Let us note that the terms which appear in this formulation are powers. The use of R e (-i ω) is not reducing, because it is redundant with the equation that would result from I m (-i ω). In practice, the VTCR procedure consists of two stages: first, one builds admissible fields of S l ad as will be described further on; then, one introduces the boundary conditions through the discretization of the variational formulation associated with a set of admissible fields chosen among the elements of S l ad . Consequently one can make a parallel with Trefftz methods [START_REF] Kita | Trefftz methods: an overview[END_REF], which also use frontier equations formulations in which fields defined on each substructure are implemented. Thus, one obtains the substructure quantity w l, n l , which must be equal to the interface quantity w l, n l . The expression of the interface/substructure l compatibility can then be summarized by the admissible space S Γ, l ad :

s l = (w l, n l , M l ) ∈ S Γ, l ad ⇔ w l, n l = f l (M l ) + g d l (7)
where g d l comes from the prescribed quantities on ∂S l , and f l is an operator coming from the condensation onto Γ.

The final problem

Once the substructure's equations have been obtained through the substructure/interface compatibility problem (7), the final step consists in solving the equations of the joints by minimizing a residual called the interface error. For example, for a perfectly rigid joint, the interface's behavior is described by the following equations:

-compatibility w 1, n 1 = w 2, n 2 -equilibrium M 1 + M 2 = 0 (8)
which is the key to the construction of the interface's L 2 residual:

E 2 Interf ace = Γ 12 w 1, n 1 -w 2, n 2 2 dL + g Γ 12 M 1 + M 2 2 dL (9)
where g must homogenizes the dimensions of the two terms. It is actually chosen in order to prevent numerical difficulties relative to the conditioning number of the problem, which are generally encountered with Trefftz methods. Then, the final problem consists in a constrained minimization:

Final Problem.

Find s l = (w l, n l , M l ) ∈ S Γ, l ad , ∀l, which minimizes:

E 2 Interf ace = Γ 12 w 1, n 1 -w 2, n 2 2 dL + g Γ 12 M 1 + M 2 2 dL (10)
One can show that this problem is equivalent to the reference problem. As will be detailed later, once the compatibility space S Γ, l ad has been discretized, this constrained minimization problem is solved using Lagrange multipliers.

Let us note that since the interfaces are treated only in a weak way and priority is given to the equations of the substructures it is important to verify that this formulation actually enable us to take into account elaborate joint models. An example will show that this is indeed the case.

Construction of the admissible fields

Let us now describe how the VTCR shape functions are defined. This stage is a fundamental and classical aspect of the method [START_REF] Ladevèze | The variational theory of complex rays for the calculation of mediumfrequency vibrations[END_REF][START_REF] Rouch | Modélisation et Calcul Des Structures Complexes En Moyennes Fréquences Par la Théorie Variationnelle Des Rayons Complexes[END_REF]. The VTCR uses a two-scale representation of subsets of S l ad . For example, in the vicinity of a point X of S l , the solution is described as the superposition of modes which can be expressed as follows:

W l (X, P ) = e i √ ωP .X U l (X, P ) (11)
The moments are deduced from the constitutive relation:

M l (X, P ) = (1 + iη l )K P S X (W l (X, P ) ) (12)
e i √ ωP .X represents the "fast" part of the field, which is defined explicitly. U (X, P ) is the "slow" part. P denotes a vector which characterizes the local vibration mode and is chosen such that the mode verifies the dynamic equations and the constitutive relation of Plate l. There are many possible options and choosing among them depends on the zone one is particularly interested in. The VTCR distinguishes the interior zone, the edge zone and the corner zone. For example, let us concentrate on the properties of n th -order complex rays for the interior zone. In this case, U l includes an n th -order polynomial in X. If the damping factor η l is small, U l can be expressed as:

U l (X, P ) = e η l 4 √ ωP .X V (X, P ) [a l (P ) ] (13)
The coefficients a l (P ) of the polynomial V (X, P ) [a l (P ) ] are unknowns; they form a column [a l (P ) ] of "generalized amplitudes". Consequently, all the unknowns are large-wavelength quantities. The damping of the wave being modeled is represented by the expression e η l 4 √ ωP .X .

These fields can be viewed as vibration modes in an infinite domain having the same mechanical properties as Substructure S l : therefore, they belong to S l ad . In the following sections, for the sake of simplicity, we will consider the particular case of a 0 th -order complex ray and omit the index l of the plate. In that case, the simplified expression of the complex ray is:

W (X, P ) = a (P ) e ( η 4 √ ωP •X) e (i √ ωP •X) = a (P ) w e (X, P ) (14)
a (P ) being the complex magnitude of the elementary complex ray w e (X, P ) . The admissibility relation

∆∆W -k 4 W = 0 on S ( 15
)
where

k 4 = 3ρω 2 (1 -ν 2 ) (1 + iη)Eh 2 requires: (P .P ) 2 = r 4 with r 4 = 3ρ(1 -ν 2 ) Eh 2 (16)
Therefore, the locus of the end of the admissible wave vector P for the interior zone is a circle C of radius r (see Fig. 3). This curve is a characteristic of the material. As one follows the circular path, every direction of the plate is taken into account. Examples of such modes can be seen in Fig. 4. The definitions of corner modes and edge modes are very similar. For the interior zone, one takes into account all admissible ray directions. The solution w sol can be expressed as:

w sol (X) = C W (X, P ) ds (17)
The integrals are calculated along Curve C. Clearly, w sol depends on the generalized amplitude [a] relative to P defined on C.

The discretized form of the VTCR

All that is needed in order to derive approximations from the VTCR is the selection of a subspace S l ad, h of S l ad . In practice, one discretizes Curve C into n elements which can be of different sizes (for example, in Fig. 5, C is discretized into eight elements of the same size.) The amplitude a (P ) is considered to be constant throughout angular sectors a h i . Then:

w h (X) = n i=1 a h i C i w e (X, P ) ds = n i=1 a h i w h i(X) = W h (X) • a h (18)
The generalized amplitude a h i is associated with the basic function w h i (X) corresponding to the integral of an elementary ray over an angular sector related to C i . Let us note that since the admissible Equations (1) are linear, the basic function w h i (X) derive from the integration in respect with the ray direction ϕ and keep on satisfying the governing equations. This point preserves the variational consistency of the VTCR approach. In practice, the integration is performed thanks to a multiscale approach : indeed, a complex rays is composed of two parts : a fast one and a slow one. The fast part of the complex ray is treated analytically while the slow part is integrated by a numerical approach.

Then, the subspace S l ad, h of the approximations deduced from S l ad is defined by:

s h l ∈ S l ad, h ⇔ w h (X) = W h (X) • a h M h (X) = Λ h (X) • a h (19)
The same procedure can be applied to corner zones and edge zones.

In order to end up with a finite-dimension problem, the interface quantities must also be discretized. Let us note that there is no obligation for the interface quantities and the substructure quantities to be discretized in the same way.

s h l ∈ S Γ, l ad, h ⇔    w h (X) = W h (X) • a h M h (X) = Λ h (X) • a h (20)
The interface/substructure compatibility problem (7) leads to the following linear problem:

(a l , a l ) ∈ S Γ, l ad, h , ∀l, ⇔ [K adm ] a a = F adm (21)
while the interface error leads to a quadratic form:

E 2 Interf ace = 1 2 t a a [K err ] a a ( 22
)
a h i w e (X, P ) a(P )w e (X, P )

Fig. 5. Admissible interior modes for a plate

Thus, the final problem (10) can be expressed as the following constrained minimization problem:

Final Discretized Problem. Find (a, a) which minimizes E 2 Interf ace = 1 2 t a a [K err ] a a under the compatibility constraints [K adm ] a a = F adm (23)
Finally, the constrained minimization problem is solved using Lagrange multipliers λ:

1 2 (K err + t K err ) -t K adm -K adm 0   a a λ   = 0 F adm (24)
Let us note that the minimization constraints do not depend on the joint parameters, which are contained in the residual: the admissibility system is built once and for all, and changes in interface behavior affect only the residual matrix K err .

Anterior studies have shown that in practice, it is not necessary to take into account order greater than zero in the slow part a h , which will be the case in the examples presented in this paper. In the rest of the presentation, for numerical reasons, we choose to discretize the interface quantities using the same space S l ad, h as for the corresponding substructure quantities, or a subspace of this space. The discretization of the interface quantities will be studied in the last part.

Effective quantities

From a mechanical point of view, the spatial distribution of the solution has no reliable physical meaning in the medium-frequency range.

Indeed, the response at a specific location and for a specific frequency is extremely sensitive to perturbations. In order to interpret the response, spatially-averaged quantities are more useful. Thus, it is necessary to extract from the solution some effective quantities defined on a domain greater than one wavelength.

On a given substructure S, let us designate the quantities of interest by Q (X) , which can represent, for example, the strain energy, the kinetic energy or the dissipation work. In the vicinity of any test point, the associated effective quantities are defined over a domain D by:

Q D = 1 D D |Q (X) | dX (25)
For instance, in the last section, the domains D considered will be composed of each plate S i of the structure.

Modeling of the interfaces

Today, the modeling of beams, plates, shells and their combinations no longer presents difficulties for the VTCR. Dedicated techniques allows us to implement more precise details like uncertainties (holes for an example) located in the interior of substructures (Ladevèze, et al. 2003). However in actual thin structures (such as cars, satellites or launchers), most of the damping, which governs the amplitude of the response to dynamic solicitations, is localized in the different connections between substructures. These joints (whether welded, riveted or bolted) modify the structural response considerably and can be used by engineers to adjust the structure's vibrational properties. As an illustration, experiments on three generic assemblies (a flat plate, a folded plate and a welded plate) clearly show that joint stiffness and damping must be taken into account. The experimental setup is not described in this paper since it is the subject of an additional one. The results are synthesized in Table I and on Fig. 6. They shows that the damping of the joints depends on the type of connexion which is used, and on the frequency of the solicitation. Consequently the modeling of the joints must be improved in order to take into account more realistic joint behavior.

Frequency (Hz)

200 -400 400 -800 800 -1600 1600 -3200 Flat plate P inj (10 The joints presented here are described by distributed interface parameters, such as mass, stiffness and damping coefficients (see Fig. 7), which enable us to model complex heterogeneous connections. These coefficients can be identified on generic assemblies to

Modeling of the joints

The local behavior can be modeled by interface relations which represent elastic and/or dissipative joint characteristics. Each joint is considered to be a mechanical entity whose own specific behavior depends on the type of connection. Many different types of joints can be modeled by this approach.

As in the case of perfect joints, there are four local boundary/interface conditions for each side of a plate: two equilibrium equations (one for the bending moments and the other for the Kirchhoff shears) and two compatibility equations (one for the displacements and the other for the slopes). For example, Equation ( 9) is the residual of perfectly rigid interface equations for bending moments: w l, n l -w m, n m = 0 continuity of the slopes M l + M m = 0 equilibrium of the bending moments (26)

The residual for Kirchhoff shears would be very similar. In order to obtain a more realistic model of the connections, let us now describe the local behavior by the following equations of a complex elastic and/or dissipative joint:

1 2 (-M l + M m ) = -k b (1 + iη b )(w l, n l -w m, n m ) constitutive relation M l + M m = -Jω 2 2 (w l, n l + w m, n m ) dynamic equation (27)
where k b and J are the line density of bending stiffness and inertia respectively, and η b models the corresponding dissipative part.

We built the residual of Equations ( 27) as:

E 2 Interf ace = Γ 12 1 2 (-M l + M m ) + k b (1 + iη b )(w l, n l -w m, n m ) 2 dL + g Γ 12 M l + M m + Jω 2 2 (w l, n l + w m, n m ) 2 dL ( 28
)
where g is chosen to balance the order of magnitude of the two terms in order to avoid numerical problems. Similarly, we can introduce new clamping conditions between the structure and the ground, which seems to be more realistic than prescribing zero slope:

M l = -k b (1 + iη b )w l, n l
bending spring with the ground (29)

This condition can be implemented thanks to the following residual:

E 2 Interf ace = Γ 0 M l + k b (1 + iη b )w l, n l 2 dL (30)
Let us note that in order to improve accuracy the interface coefficients, such as stiffness, mass and damping, can be made to vary along the interface. For example, a riveted joint can be modeled by a stiffness coefficient which is infinite and zero alternatively, as shown in Fig. 8 (where the joint's stiffness is denoted k). In order to implement such parameters numerically, it is necessary to discretize them. Considering that in actual structures most of the joints have periodic characteristics, we chose to use a Fourier decomposition of the interface coefficients, such as:

k (X) = +∞ n=-∞ k n e i n λ•X +N n=-N k n e i n λ•X X ∈ Γ (31)
with λ being the period of the pattern describing the junction. Fig. 8 and Fig. 9 show an example of such a calculation in which the joint's bending stiffness and translation stiffness have both been taken into account. The VTCR calculations Fig. 8. The heterogeneous interface model and its calculation required 260 degrees of freedom (128 interior modes and 33 edge modes per edge). The results confirm that the boundary conditions are accurately respected. In this case, modeling the joint with constant parameters along the interface could not have given an acceptable result. In order to validate this VTCR solution, we compared it with the solution obtained with the traditional FEM code NASTRAN. The FEM calculation was conducted using at least 10 elements per wavelength, according to the so-called "rule of the thumb" (see [START_REF] Ihlenburg | Finite Element Solution of the Helmholtz Equation with High Wave Number. Part I: The h-Version of the FEM[END_REF], [START_REF] Ihlenburg | Finite Element Solution of the Helmholtz Equation with High Wave Number. Part 2: The h-p-Version of the FEM[END_REF]). 82, 863 degrees of freedom were required. The maximum displacements obtained were 5.26 10 -4 m with the VTCR and 4.96 10 -4 m with NASTRAN, which amounts to a difference of about 6 %.

Numerical example

In this section, we present the influence of the joint's bending stiffness and damping on the response of an assembly of two plates connected by a complex joint (see Fig. 10). Here, the interface parameters were considered to be constant along the interface. To make comparisons possible, we used the effective displacements of each plate as defined in 2.5. The domains D considered in this definition are composed of each plate S i of the structure. These quantities were averaged over [100Hz,200Hz]. Let us note that for the calculations with the VTCR only 200 degrees of freedom were used: 68 d.o.f. for a plate (32 interior modes and 9 edge modes per edge), and 32 other d.o.f. for the interface quantities relative to the plate. The minimization contraints add 68 Lagrange multipliers per plate (as many as the number of test functions of each plate). Consequently, the final size of the problems is 336 degrees of freedom, which seems to be very reasonable.

We can observe on Fig. 11 that a stiffness per unit length greater than 10 3 N.m.rd -1 is equivalent to a perfectly rigid joint. Similarly, for a joint stiffness equal to 10 2 N.m.rd -1 , damping is significant if it is greater than 10 -1 . This value can seem very high, but it should not be compared directly to the damping of the plate (equal to 0.01) because what matters is the contribution of the joint's dissipation to the total energy dissipation of the structure.

Example

Let us compare the results obtained with the classical VTCR formulation, used as the reference, with those obtained by the new formulation using the minimization of the interface's residual. Here the discretization on each plate required 68 degrees of freedom for the classical VTCR formulation (32 interior basis functions per plate, 9 edge modes per edge), and the new formulation introduces for each plate 32 interface degrees of freedom and 68 Lagrange multipliers respectively. Fig. 12 shows the solutions obtained by the two methods for the two-plate assembly with a perfectly rigid joint and a 220 Hz solicitation. The solutions are very similar, as can be seen from the effective displacements (Table II The next calculations concern the same assembly of two plates, solicited at 100 Hz, but with a complex joint and different values of the joint's stiffness. The discretization requires 68 degrees of freedom per plate (32 interior modes and 9 edge modes per edge) for each substructures. Different discretizations were used for the interface quantities: discretizations (a), (b) and (c) consisted of the border complex rays (9 edge modes per edge), the interior complex rays (32 interior modes per plate), and both the border and the interior complex rays (32 interior modes and 9 edge modes per edge) respectively. (d) is the reference calculated with the classical global formulation involving 68 degrees of freedom per plate (32 interior modes and 9 edge modes per edge).

Fig. 13 shows that the results are exactly the same and that a coarse discretization of the interface is sufficient, provided that it contains the interface's complex rays relative to the joint being considered.

Conclusion

The Variational Theory of Complex Rays is a predictive approach which is very well-suited to medium-frequency vibration analysis of complex structures. Its main features are the use of shape functions with a strong mechanical meaning and a variational formulation which results in a very low cost, independent of the frequency. The variational formulation is associated with the boundary conditions, which is very useful in order to take into account different complex connections between substructures. Since the modeling of actual heterogeneous joints is very important in order to obtain an accurate response in mediumfrequency vibrations with this type of structure, we introduced more sophisticated interface equations through mass-stiffness-damping joint behavior. We presented a new formulation of the theory which enables one to analyze the same structure with different connections more conveniently. This new formulation is based on the minimization of a L 2 residual of interface equations while the admissibility of the well-modeled substructures is enforced one and for all by minimization constraints. This formulation should not meet any theoretical difficulty to treat more complex structures like assemblies of beams, plates, or shells. In future works, we will take advantage of this formulation to identify the parameters of the joints using experimental results.

Fig. 3 .

 3 Fig. 3. Admissible P for interior modes

Fig. 6 .

 6 Fig. 6. Global damping values for several frequency ranges

Fig. 7 .

 7 Fig. 7. Generic assembly modeling

Fig. 9 .

 9 Fig. 9. Heterogeneous interface (245 Hz): (a) VTCR solution, (b) NASTRAN solution

Fig. 10 .

 10 Fig. 10. Assembly of two plates with a complex constant joint

Fig. 12 .

 12 Fig. 12. Solutions obtained with the two formulations

Fig. 13 .

 13 Fig. 13. Comparison between the classical VTCR and the substructured VTCR

Table II .

 II ). Validation of the substructured version of the VTCR

		Effective quantities
		Classical formulation Substructured formulation
	Plate 1	16.333 10 5	16.389 10 5
	Plate 2	12.117 10 5	12.061 10 5

O. Dorival, P. Rouch, O. Allix

O. Dorival, P. Rouch, O. Allix


Olivier Dorival
email: dorival@lmt.ens-cachan.fr

Philippe Rouch
email: rouch@lmt.ens-cachan.fr

Olivier Allix
email: allix@lmt.ens-cachan.fr

A method for updating joint parameters in medium-frequency vibrations

Keywords: medium frequency, joints, damping, updating, VTCR, inverse problem, domain decomposition methods moyennes fréquences, liaisons, amortissement, recalage, TVRC, problème inverse, décomposition de domaine

Joints between substructures play a significant role in the vibrational behavior of complex structures because they govern energy flow and most of the dissipative phenomena. In order to identify joint models, this paper proposes a robust updating method which was initially based on studies of the error in constitutive relation in relation to finite element model updating. Here, it is redesigned in order to focus on joint models in medium-frequency problems. In order to do that, we use an alternative numerical approach called the Variational Theory of Complex Rays (VTCR). After introducing the new formulation, the paper analyzes the effectiveness of the approach in identifying a joint's stiffness and damping.

RÉSUMÉ. Les liaisons entre sous-structures jouent un rôle crucial dans la réponse vibratoire d'une structure complexe, en régissant les transferts énergétiques ainsi qu'une grande partie de la dissipation. Pour identifier ces modèles de liaison, ce travail s'inspire d'une méthode de recalage robuste, initialement basée sur les travaux sur l'erreur en relation de comportement pour la correction de modèles éléments finis. La formulation est revisitée pour se concentrer sur le problème des jonctions en moyennes fréquences. Das ce but une approche numérique alternative aux éléments finis, la théorie variationnelle des rayons complexes (TVRC) est utilisée. Après avoir présenté la nouvelle formulation, nous étudions l'efficacité de cette approche pour l'identification des paramètres de raideur et d'amortissement d'une liaison.

Introduction

The treatment of dissipation is of primary importance in seeking the numerical response of vibrating structures. In many real systems, one can decompose the structure into simple linear viscoelastic substructures and complex, coarsely modeled interfaces. The behavior of joints, which usually involves complex dissipative phenomena (micro/macro frictional contact, microshocks, heating, ...), must be accounted for because in vibration problems energy losses affect the magnitude of the structural response directly. In addition, the stiffness of joints can significantly influence how the injected power is distributed among the different substructures. Regarding the low-frequency range, updating methods for joints have already been the subject of much research [START_REF] Mottershead | Model updating in structural dynamics: a survey[END_REF][START_REF] Deraemaeker | Sur la maîtrise des modèles en dynamique des structures à partir de résultats d'essais[END_REF]; for high frequencies, the identification of coupling loss factors [START_REF] Maxit | Estimation of SEA coupling loss factors using a dual formulation and FEM modal information, Part I: theory[END_REF] is an unavoidable step of methods derived from Statistical Energy Analysis (SEA) [START_REF] Lyon | Statistical Energy Analysis of dynamical systems: theory and applications[END_REF]. Between these two extremes, medium frequencies present major difficulties because of their high modal density and their high sensitivity to boundary conditions. This paper presents a method for updating joint parameters -in our case, joint stiffness and damping -with special emphasis on the robustness of this technique.

Different methods are available in order to deal with various types of inverse physical problems [START_REF] Bonnet | Inverse problems in elasticity[END_REF]. Our method uses the modified error in constitutive relation approach [START_REF] Ladevèze | Application of a posteriori error estimation for structural model updating[END_REF], which is derived from the well-known concept of error in constitutive relation for the quantification of finite element discretization errors [START_REF] Ladevèze | Error estimate procedure in the finite element method and application[END_REF]. Starting with a complete set of available information coming from the numerical model as well as experimental measurements, the method consists in dividing the corresponding equations into two groups according to their reliability. The reliable equations (in relation to both the numerical model and the experimental model) are to be satisfied exactly, while the unreliable equations can be solved approximately through the minimization of a cost function. In the case of a joint modeling problem, the equations of the continuous model are divided as follows: the substructures are considered to be reliable components, while the joint equations -including the model parameters which need to be identified -are unreliable parts of the structure. Regarding the experimental information, the locations and directions of the sensors and the prescribed angular frequency and excitation force are assumed to be reliable. The quantities measured by the sensors, which are unavoidably noisy, are considered unreliable. Table 1 summarizes the characteristics of these different types of information. The solution of the unreliable equations is achieved by minimizing a cost function composed of two terms: a modeling error which quantifies the satisfaction of the joint equations, and an experimental error which represents the discrepancy between the numerical solution and the experimental measurements.

Then, a numerical scheme must be applied in order to solve the problem. In dealing with medium-frequency problems, classical methods encounter serious difficulties. Finite element calculations require the use of a very refined mesh in order to avoid what one calls pollution errors due to the small-wavelength phenomena involved [START_REF] Ihlenburg | Finite element solution of the Helmholtz equation with high wave number. Part II: The h-p-version of the FEM[END_REF]; this leads to sharp increases in computational costs.

High-frequency methods must be modified in order to release some of the SEA assumptions which are no longer valid in the medium-frequency regime. Among the various possible methods [START_REF] Farhat | A discontinuous Galerkin method with Lagrange multipliers for the solution of Helmholtz problems in the mid-frequency regime[END_REF][START_REF] Strouboulis | The Generalized Finite Element Method for Helmholtz equation: theory, computations and open problems[END_REF], we chose to use the Variational Theory of Complex Rays (VTCR) [START_REF] Ladevèze | Analysis of medium-frequency vibrations in a frequency range[END_REF], which has been found to be easy to adapt to the separation of the equations discussed previously. This Trefftz method consists in building over each substructure two-scale shape functions which satisfy the dynamic equilibrium and the constitutive law of the substructure exactly. (These equations are considered to be reliable.) The boundary conditions of each subdomain are applied weakly through a variational equation which can be viewed as a balance between virtual powers at the boundaries. In order to distinguish joint equations from other boundary conditions, we use a substructured formulation of the VTCR which was presented in [START_REF] Dorival | A substructured version of the variational theory of complex rays dedicated to the calculation of assemblies with dissipative joints in the mediumfrequency range[END_REF]. Further details of the formulation and its discretization can be found in [START_REF] Dorival | A substructured Trefftz method for updating joint models in the medium-frequency range[END_REF]. This paper is organized as follows: Section 2 briefly describes the formulation. Section 3 presents the identification of the joint's stiffness and damping parameters for a simple simulated structure, which illustrates the robustness of the method. -Measured amplitudes Unreliable -constitutive law of the joint (with noise) information -dynamic equilibrium of the joint 2. Formulation of the inverse problem using a substructured Trefftz method

A substructured Trefftz method for the equations of the continuous model

This first section deals with the equations of the continuous model alone. The experimental equations will be introduced in the next section in order to produce an inverse problem. The approach is based on a substructured version of the VTCR which was detailed in [START_REF] Dorival | A substructured version of the variational theory of complex rays dedicated to the calculation of assemblies with dissipative joints in the mediumfrequency range[END_REF]. For the sake of simplicity, let us consider two substructures Ω i, i=1, 2 and denote the displacement field and the stress tensor u i and σ i respectively.

The substructures Ω i are distinct from the joint Γ to be identified (see Figure 1). The joint is considered to be a real substructure with its own equations and interface unknowns denoted Ŵ i and F i (which are distinct from the unknowns of the adjacent substructures). Interface forces F i are applied to each substructure. The equations of the substructures are considered to be reliable and satisfied exactly. The numerical treatment of each substructure Ω i involves a classical formulation of the VTCR: the solution within each substructure is sought as a linear combination of complex rays (see Figure 2) which satisfies the local dynamic equilibrium and constitutive law exactly. Boundary conditions are applied in variational form. In practice, the strong physical meaning of the shape functions enables one to use only a very small number of unknowns -about one hundred per substructure. It is important to observe that the numerical cost does not increase with the frequency being considered. Additional details on the VTCR can be found in [START_REF] Ladevèze | Analysis of medium-frequency vibrations in a frequency range[END_REF]. The problems in each substructure Ω i lead to admissibility constraints associated with admissibility spaces Joint equations with joint unknowns Ŵ i and F i must be satisfied. In our case, we assume that line spring-mass-dampers are distributed along the interface Γ. The equations are satisfied by minimizing a residual, called the modeling error, which is constructed as follows:

S i ad . Γ F 2d u 1d Ω 1 Ω 2 ∂ F Ω 2 ∂ u Ω 1 ∂ F Ω 1 ∂ u Ω 2 u 2d F 1d W 2 W 1 (u 1 , σ 1 ) (u 2 , σ 2 ) F 2 F1
E 2 mdl = 1 NN Γ F 1 + F 2 + m ω 2 2 ( Ŵ 1 + Ŵ 2 ) 2 dl + 1 NW Γ Ŵ 1 -Ŵ 2 + k 2 ( F 1 -F 2 ) 2 dl [1]
The line mass m and line stiffness k can vary along the joint and dissipative behavior can be represented as a complex part of m and k. Then, the joint unknowns are discretized in the same way as the matching substructure unknowns in order to obtain a conforming discretization.

The continuous problem consists in minimizing the modeling error E 2 mdl under the admissibility constraints S i ad . This is achieved by solving the corresponding lagrangian problem.

Experimental information and the inverse problem

In order to produce an inverse problem, experiments are performed on the actual structure. Let us designate measured quantities as •. When dealing with mediumfrequency vibrations, the hypersensitivity of the structure requires the use of nonclassical measurement techniques. Indeed, intrusive techniques modify the structure's mass and stiffness locally, and pointwise measurements provide only very local information. Laser vibrometry, which is widely used in high-frequency measurements, is an appealing technique to overcome these difficulties because it is nonintrusive and can provide field measurements capable of capturing the very local phenomena which occur in this frequency range.

Given a set of measurement points X m , let us define the following measurement error in order to quantify the distance between the simulation and the experimental results:

E 2 mes = 1 N mes m u (X m ) -u (X m ) 2 [2]
The summation of the measurement error and the modeling error leads to the modified error:

E 2 mod = (1 -r) E 2 mdl + r E 2 mes [3]
where r is used to balance the two terms according to one's relative confidence in the model and in the experiments. If r = 0 is chosen, the modified error only accounts for the interface model. On the contrary if r = 1, the interface equations have no weight and the solution tries to satisfy the experimental results at best. For r chosen between 0 and 1, one obtains a solution which is a compromise between the measured data and the joint model chosen by minimizing the modified error under admissibility constraints S i ad , i = 1, 2 . A typical value of 0.8 will be used in Section 3. Since the modified error depends on k joint parameters p k , the inverse problem consists in finding the optimum joint parameters p opt k :

(p opt k ) = arg min (p k ) E 2 mod (p k ) [4]
The corresponding nonlinear optimization problem can be solved through classical Newton-like methods. For the sake of simplicity, the shear behavior of the joint is assumed to be perfectly rigid and, therefore, the rotational behavior of the joint alone needs to be identified. In our case, the joint's rotational stiffness of the "actual" structure was k 0 = 10 1.5 N.m and its mass was neglected. Viscous joint damping was assumed with two dissipation levels: a slightly dissipative joint with η = 0.1 and a highly dissipative joint with η = 10. At this frequency, the first case dissipated 9 % of the total energy and the second case 68 %. The rotational joint equations are:

n 1 M1 n 1 + n 2 M2 n 2 = 0 moment equilibrium n M n = -k 0 (1 + iη) ū,n constitutive law [5]
where Mi and ûi ,n respectively denote the bending moment and the normal rotation over Γ. The overbar notation ū,n (respectively M) represents the average value of normal rotations (resp. bending moments) over Γ. In order to simulate measurement noise, the experimental data were perturbed by ± 20 % uniform random noise.

In the identification process, a classical value of 0.8 was assigned to the weighting parameter r. The inverse problem was discretized using 400 complex rays -64 interior complex rays and 4 x 9 edge rays for each plate, and the same for the interface DOFs. In the following section, we investigate the possibility of identifying the joint's stiffness and damping separately. Indeed, a real structure can require the identification of a large number of joints, and the optimization problem [4] can become very costly. Therefore, separate identification can be useful. In addition, for identifiability reasons [START_REF] Wang | Experimental identification of mechanical joint parameters[END_REF], it is advisable to seek the stiffness first, then the damping coefficient. The effectiveness of this approach is studied in the next section.

Identification of joint stiffness

(N.m) Modified error E E mod -exact η E mod -η +50% E mod -η -50% E mod -η +300% E mod -η -95 % exact k 0 2 2 2 2 2 2 mod Noise 20% + - r = 0.8 Figure 4.
The negligible influence of an erroneous η on the minimum of the cost function E 2 mod (k0) for a slightly dissipative joint. The vertical dotted line represents the exact parameter to be identified (where the minimum should be located)

Figure 4 shows the evolution of the modified error when identifying the stiffness of a slightly dissipative joint (η = 0.1). The cost function is sufficiently convex to be minimized without difficulty both for exact and highly erroneous damping coefficients. The errors on η have no impact: the identified stiffness is very close to the exact stiffness represented by a vertical dotted line. This is easily explained by the fact that for this slightly dissipative joint (in comparison with the substructures' dissipation) the joint's damping coefficient does not play a significant role. Let us note that the minimum of the cost function is due to the measurement noise, which has no influence on the location of the minimum because of the design of the cost function and the field measurements.

Figure 5 shows the same results for a highly dissipative joint (η = 10). With an exact damping coefficient (continuous line), the joint's stiffness is correctly identified. This is not true with an erroneous damping coefficient, in which case the identified stiffness differs from the exact value (vertical dotted line). We conclude that for a highly dissipative joint the damping parameter η plays a significant role and must be modeled accurately in order to recover the correct stiffness. Modified error mod (k0) for a highly dissipative joint. The vertical dotted line represents the exact parameter to be identified (where the minimum should be located)

E E mod -exact η E mod -η +50% E mod -η -50% E mod -η +300% E mod -η -95 % exact k 0

Identification of the joint damping

In this section, we address the identification of the joint's damping with particular emphasis on the influence of an erroneous joint stiffness. To understand this better, the measurement noise was removed. The evolution of the modified error as a function of the damping parameter η in the case the highly dissipative joint is shown in Figure 6. The minimization of the cost function leads to the exact damping parameter, even if the stiffness is highly erroneous. The explanation is quite simple: in the case of a highly dissipative joint (in comparison with the substructures' dissipation), the response is influenced mainly by the damping coefficient, which leads to a very robust identification of η.

E mod -exact k 0 E mod -k 0 +10% E mod -k 0 -10% E mod -k 0 +50% E mod -k 0 -50% exact η
E mod -exact k 0 E mod -k 0 +10% E mod -k 0 -10% E mod -k 0 +50% E mod -k 0 -50% exact η 2 2 2 2 2 2 Unnoised r = 0.8 mod Figure 7.
The strong influence of an erroneous k 0 on the minimum of the cost function E 2 mod (η) for a slightly dissipative joint. The vertical dotted line represents the exact parameter to be identified (where the minimum should be located)

In the case of a slightly dissipative joint (see Figure 7), the evolution of the cost function is less affected by the damping parameter η. If the stiffness is taken equal to the exact value (continuous line), the damping parameter identified at the minimum of the cost function is exact. Unfortunately, the existence of erroneous joint stiffnesses perturbs this minimum, or even causes it to vanish, which makes the updating process even trickier. In this case, the influence of the damping parameter is not significant enough because the dissipation in the joint is negligible compared to the dissipation in the substructures. Consequently, classical difficulties in the identification of η [START_REF] Wang | Experimental identification of mechanical joint parameters[END_REF] seem unavoidable.

Concluding remarks

In this paper, we presented an updating method for joint models. This technique uses a suitable separation of all information in order to circumvent the difficulties related to medium-frequency vibrations as well as emphasize the behavior of the joints. The numerical approach for discretizing the inverse problem follows exactly the same lines. The robustness of the method was briefly illustrated using a simulated structure. We found that the identification of a joint's stiffness and damping by a step-by-step process cannot be successful if the joint's dissipation is high. In this case, damping should be identified first. This conclusion is a consequence of the sensitivity of the response to the two parameters. At this stage, the experience of the engineer is of great importance for industrial applications. Although the numerical framework developed here turns out to be effective and efficient, further studies should focus on the derivation of a physically consistent experimental part (in the cost function) suitable for medium-frequency phenomena.

Figure 1 .Figure 2 .

 12 Figure

Figure 3 .

 3 Figure 3. Simulated experimental structure consisting of two plates connected by a joint

Figure 5 .

 5 Figure 5. The strong influence of an erroneous η on the minimum of the cost function E 2mod (k0) for a highly dissipative joint. The vertical dotted line represents the exact parameter to be identified (where the minimum should be located)

Figure 6 .

 6 Figure6. The negligible influence of an erroneous k 0 on the minimum of the cost function E 2 mod (η) for a highly dissipative joint. The vertical dotted line represents the exact parameter to be identified (where the minimum should be located)

Table 1 .

 1 The proposed method for carrying out joint updating

		Continuous model	Experimental model
		-Geometry	-Angular frequency
	Reliable	-Substructure equations:	-Sensor locations
	Information	-local dynamic equilibrium	-Sensor directions
		-constitutive law	-Prescribed forces
		-Reliable boundary conditions
		-Joint equations: