Nobutaka Ito
email: ito.nobutaka@lab.ntt.co.jp

Emmanuel Vincent

Tomohiro Nakatani

Nobutaka Ono

Shoko Araki

Shigeki Sagayama

Blind suppression of nonstationary diffuse noise based on spatial covariance matrix decomposition

We propose methods for blind suppression of nonstationary diffuse noise based on decomposition of the observed spatial covariance matrix into signal and noise parts. In modeling noise to regularize the ill-posed decomposition problem, we exploit spatial invariance (isotropy) instead of temporal invariance (stationarity). The isotropy assumption is that the spatial cross-spectrum of noise is dependent on the distance between microphones and independent of the direction between them. We propose methods for spatial covariance matrix decomposition based on least squares and maximum likelihood estimation. The methods are validated on real-world recordings.

Introduction

Noise suppression is the task of estimating a desired signal from its noisy observations by microphones. The difficulty of the noise suppression problem depends on the temporal and the spatial properties of noise. We can distinguish four categories of noise based on two independent axes: stationary or nonstationary noise, and point-source or diffuse noise. Here, diffuse noise is defined as noise from numerous directions caused by many point sources (e.g., many interfering speakers) or by a continuous source (e.g., the vibrating body of a vehicle).

Most of the conventional approaches to noise suppression have assumed either stationary or point-source noise. The spectral approach [START_REF] Boll | Suppression of acoustic noise in speech using spectral subtraction[END_REF][START_REF] Martin | Spectral subtraction based on minimum statistics[END_REF][START_REF] Ephraim | Speech enhancement using a minimum-mean square error short-time spectral amplitude estimator[END_REF] can suppress stationary noise. Under the assumption of stationary noise, the signal power spectrum and/or the noise power spectrum are estimated and used to design spectral filters. On the other hand, the spatial approach [START_REF] Dudgeon | Array Signal Processing: Concepts and Techniques[END_REF][START_REF] Brandstein | Microphone Arrays: Signal Processing Techniques and Applications[END_REF] can suppress point-source noise by controlling the directivity, but cannot suppress diffuse noise sufficiently. This is because the number of spatial nulls is limited to the number of microphones minus one. These conventional approaches do not cover the remaining category: nonstationary, diffuse noise, which is omnipresent in the real world. Indeed, such noise is encountered in many environments, such as stations, airports, vehicles, factories, cafeterias, bars, streets, etc. This has significantly limited the application area of noise suppression techniques to the real world.

In principle, even nonstationary diffuse noise can be suppressed with a spectral filter, provided that the signal is sufficiently sparse in the time-frequency domain. Indeed, if we attenuate energy at the time-frequency points where the signal is inactive, the signal-to-noise ratio is expected to improve significantly. The question is how we can design such a filter without assuming noise stationarity. Since this is a highly ill-posed inverse problem, we need to restrict the search space by modeling the signal and noise appropriately.

To this end, we model noise based on spatial invariance (isotropy) instead of temporal invariance (stationarity). The isotropy assumption is that the spatial cross-spectrum of noise is dependent on the distance between microphones and independent of the direction between them. We propose methods for spatial covariance matrix decomposition based on least squares (LS) and maximum likelihood (ML) estimation. The estimated spatial covariance matrices of the signal and noise are used to design a time-varying multichannel Wiener filter for diffuse noise suppression. We discuss several variants of the isotropic noise model that correspond to additional assumptions about the noise properties or the array geometry. We also propose a general linear algebraic framework for treating these models in a unified manner. This allows us to derive general algorithms each to all noise models, instead of different algorithms for each noise model.

Both the LS and the ML methods have pros and cons. First, as shown later, ML is equivalent to the minimization of the Itakura-Saito divergence, which has proven to be effective for audio signal processing [START_REF] Itakura | Analysis synthesis telephony based on the maximum likelihood method[END_REF]. Second, based on a probabilistic generative model, the ML method can be combined 1) with other signal processing techniques based on a generative model [START_REF] Duong | Underdetermined reverberant audio source separation using a full-rank spatial covariance model[END_REF][START_REF] Nakatani | Speech dereverberation based on variance-normalized delayed linear prediction[END_REF][START_REF] Vincent | From blind to guided audio source separation[END_REF] to deal with more general environments, and 2) with prior distributions to improve performance using our prior knowledge. Third, in terms of the computational time, the ML method is more expensive than the LS method, because the former employs matrix inversion at each iteration.

The rest of this paper is structured as follows. In Section 2, we describe our observation model and the multichannel Wiener filter for diffuse noise suppression. In Section 3, we describe our framework for modeling the signal and the diffuse noise. Section 4 proposes methods for spatial covariance matrix decomposition based on least squares estimation. Section 5 proposes a method for spatial covariance matrix decomposition based on maximum likelihood estimation. We evaluate the proposed methods experimentally in Section 6, and conclude in Section 7.

Background

In this section, we first describe our observation model. We then derive the time-varying multichannel Wiener filter for diffuse noise suppression, and point out that the design of the filter boils down to the estimation of the spatial covariance matrices of the signal and diffuse noise.

Observation model

In this paper, we focus on blind enhancement of a single desired signal in the presence of diffuse noise. We assume that the source location is stationary. Extension to the case of multiple and/or moving sources is part of future work.

In this paper, unless we note otherwise, we represent signals in the time-frequency domain (e.g., in the short-time Fourier transform domain). We denote the number of frames by T , and the frame index by t ∈ {1, . . . , T }. We omit the frequency bin index for brevity. This should not cause confusion, since each frequency bin is processed independently in this paper. The observed signal can be modeled by

y t = x t + v t , ( 1
)
where the variables are defined as follows:

-

y t ∈ C M : M -channel observed signal, -x t ∈ C M : M -channel target signal, -v t ∈ C M : M -channel diffuse noise.
The problem of diffuse noise suppression considered in this paper is formally defined as the estimation of X {x t } 1≤t≤T from Y {y t } 1≤t≤T . Here, the notation {y t } 1≤t≤T stands for {y t |1 ≤ t ≤ T }, for example.

For simplicity, we make the following assumptions on x t and v t :

-Temporal independence: {x t } 1≤t≤T is an independent series. That is, for different t, u, x t and x u are independent. {v t } 1≤t≤T is also an independent series. -Mutual independence: {x t } 1≤t≤T and {v t } 1≤t≤T are mutually independent. That is, for any t, u, x t and v u are independent. -Gaussianity: x t and v t are zero-mean complex-valued Gaussian variables with covariance matrices

Φ x t E[x t x H t ] and Φ v t E[v t v H t ].
Here, the probability density function of the complexvalued Gaussian distribution with mean µ ∈ C M and covariance matrix Σ ∈ C M ×M is given by

N C (x; µ, Σ) 1 π M det Σ exp [ -(x -µ) H Σ -1 (x -µ) ] .
(

) 2
Modeling temporal correlation due to reverberation and more complex distributions are part of future work.

Time-variant multichannel Wiener filter for diffuse noise suppression

Here we review the time-varying multichannel Wiener filter for diffuse noise suppression. The filter is a timevarying spatiotemporal filter that is optimal in the sense of linear minimum mean square error (LMMSE). Designed properly, it can suppress diffuse noise effectively, by attenuating energy at the time-frequency points at which the signal is inactive. The LMMSE estimator is defined as the linear estimator of the form

xt = W H t y t (3)
that minimizes the mean square error

E[∥x t -xt ∥ 2 ]. (4)
In this paper, we denote the Frobenius/Euclidean norm of a matrix/vector by ∥•∥, and the Hermitian transpose by the superscript H . Partial differentiation of (4) w.r.t. W * t reveals that the optimum estimator is given by [START_REF] Duong | Underdetermined reverberant audio source separation using a full-rank spatial covariance model[END_REF][START_REF] Simmer | Post-filtering techniques[END_REF][START_REF] Doclo | GSVD-based optimal filtering for single and multimicrophone speech enhancement[END_REF]

xt = Φ x t (Φ x t + Φ v t ) -1 y t , (5)
where * denotes complex conjugation. The time-varying multichannel Wiener filter ( 5) is also optimal in the senses of maximum a posteriori (MAP) and minimum mean square error (MMSE), when the signal and noise are modeled as zero-mean Gaussian random variables (see (37)).

Spatial covariance matrix decomposition

To design the time-varying multichannel Wiener filter [START_REF] Brandstein | Microphone Arrays: Signal Processing Techniques and Applications[END_REF], we need to know the spatial covariance matrices of the signal and the noise, Φ x t and Φ v t . However, we are not given these matrices, but only the spatial covariance matrix of the observed noisy data, Φ y t E[y t y H t ]. Therefore, to suppress diffuse noise effectively, it is crucial to estimate Φ x t and Φ v t accurately from Φ y t . Under the assumption that the signal and the noise are mutually uncorrelated, these matrices are related as

Φ y t = Φ x t + Φ v t . ( 6
)
Hence, we call the estimation of Φ x t and Φ v t spatial covariance matrix decomposition, which is the main focus of the paper.

Spatial covariance matrix decomposition is a highly ill-posed inverse problem. Indeed, we need to estimate both Φ x t and Φ v t given Φ y t only. Without any constraints on Φ x t and Φ v t , there would be infinitely many decompositions. Therefore, to obtain a reliable decomposition, it is important to restrict the search space by modeling Φ x t and Φ v t appropriately. Regarding the target signal, the assumption that it is emitted from a stationary (i.e., not moving) point source implies that its spatial characteristics do not change over time. Therefore, we can consider the following model of Φ x t [START_REF] Duong | Underdetermined reverberant audio source separation using a full-rank spatial covariance model[END_REF]:

Φ x t = ϕ x t B x . (7)
Here, ϕ x t ∈ R is a time-varying parameter representing the power spectrum, and B x ∈ C M ×M is a constant parameter called the coherence matrix representing the spatial characteristics. Especially, for low reverberation, B x can be approximated as a rank-one matrix

B x = hh H , ( 8
)
where h is called the steering vector. More generally, we can model the eigenvalues of B x to be sparse.

Matrix linear subspace for unified treatment of diffuse noise models

In contrast to the target signal, diffuse noise originates from many sources or from a continuous source, and therefore its spatial characteristics vary over time. For the case of many sources, for example, even if they are not moving, the spatial characteristics of the whole noise vary over time, since different sources are active at different time-frequency points. Furthermore, diffuse noise is also spectrally nonstationary in general (e.g., consider many concurrent speakers at a cocktail party).

To this issue, in [START_REF] Ito | Robust Microphone Array Signal Processing against Diffuse Noise[END_REF][START_REF] Ito | General algorithms for estimating spectrogram and transfer functions of target signal for blind suppression of diffuse noise[END_REF], we have proposed diffuse noise modeling based on isotropy. The isotropy assumption is that the spatial cross-spectrum of noise is dependent on the distance between microphones and independent of the direction between them. Under this assumption, we can show that Φ v t belongs to a low-dimensional subspace V of the vector space H over R spanned by M ×M Hermitian matrices [START_REF] Ito | Robust Microphone Array Signal Processing against Diffuse Noise[END_REF][START_REF] Ito | General algorithms for estimating spectrogram and transfer functions of target signal for blind suppression of diffuse noise[END_REF]. We call V a matrix linear subspace because it is a subspace of the space of matrices, H.

There are several choices of V depending on additional assumptions regarding the array geometry or the noise field. Firstly, the most general choice with no additional assumptions is the real-valued noise covariance model. The isotropy assumption implies that Φ v t is symmetric, and covariance matrices are Hermitian by definition. Therefore, it follows that Φ v t is real-valued:

V real = { Φ v t ∈ H Φ v t ∈ R M ×M } . ( 9
)
Secondly, for spatial class of array geometries called crystal arrays, the blind noise decorrelation (BND) model can be applied [START_REF] Ito | Robust Microphone Array Signal Processing against Diffuse Noise[END_REF][START_REF] Ito | Diffuse noise suppression using crystal-shaped microphone arrays[END_REF]:

V BND = { Φ v t ∈ H U H Φ v t U is diagonal } . ( 10
)
Here, U is a known constant unitary matrix. Thirdly, if we assume that diffuse noise is uncorrelated between two points in the limit of a large distance between them, the following spatially uncorrelated noise model holds for large arrays [START_REF] Zelinski | A microphone array with adaptive postfiltering for noise reduction in reverberant rooms[END_REF]:

V uncor = { Φ v t ∈ H Φ v t is diagonal } . ( 11
)
Finally, assuming that the microphone coordinates are known and that noise planewaves with the same power spectrum arrive from all directions, the following fixed noise coherence model holds [START_REF] Mccowan | Microphone array postfilter based on noise field coherence[END_REF]:

V coh = { ϕ v t B v ϕ v t ∈ R } . ( 12
)
Here, B v is called a noise coherence matrix, and its (m, n)th entry is given by

b v mn = sinc ( 2πf L mn c ) . ( 13
)
Here, sinc(•) denotes the sine cardinal function, L mn the distance between the mth and the nth microphones, and c the sound velocity. Each of these noise models was applied to non-blind diffuse noise suppression with a known target steering vector [START_REF] Ito | Diffuse noise suppression using crystal-shaped microphone arrays[END_REF][START_REF] Zelinski | A microphone array with adaptive postfiltering for noise reduction in reverberant rooms[END_REF][START_REF] Mccowan | Microphone array postfilter based on noise field coherence[END_REF][START_REF] Ito | Designing the Wiener post-filter for diffuse noise suppression using imaginary parts of inter-channel cross-spectra[END_REF], and to noise-robust multiple source localization [START_REF] Ito | Crystal-MUSIC: Accurate localization of multiple sources in diffuse noise environments using crystal-shaped microphone arrays[END_REF]. In [START_REF] Ito | Diffuse noise suppression using crystal-shaped microphone arrays[END_REF], we proposed the first method for blind diffuse noise suppression, which was limited to the blind noise decorrelation model. In this paper, we propose general blind diffuse noise suppression methods applicable to all the above noise models.

Spatial covariance matrix decomposition based on least squares estimation

Cost function: square error

In this section, we propose two methods for spatial covariance matrix decomposition based on the minimization of the following Euclidean error subject to Φ v t ∈ V:

J LS (Θ) ∑ t ∥Φ y t -(Φ x t + Φ v t )∥ 2 (14) = ∑ t ∥Φ y t -(ϕ x t B x + Φ v t )∥ 2 . ( 15
)
Θ denotes the parameter set given by

Θ {{Φ x t } 1≤t≤T , {Φ v t } 1≤t≤T } (16) = {{ϕ x t } 1≤t≤T , B x , {Φ v t } 1≤t≤T }. ( 17
)
The minimizers of [START_REF] Zelinski | A microphone array with adaptive postfiltering for noise reduction in reverberant rooms[END_REF] are not unique, and so an additional constraint is necessary. To see this, note that ( 15) is decomposed into V and V ⊥ components as follows:

J LS (Θ) = T ∑ t=1 ∥P[Φ y t ] -ϕ x t P[B x ] -Φ v t ∥ 2 (18)
+ T ∑ t=1 ∥P ⊥ [Φ y t ] -ϕ x t P ⊥ [B x ]∥ 2 .
Here, V ⊥ denotes the orthogonal complement of V, and P and P ⊥ the orthogonal projection operators onto V and V ⊥ , respectively. We can eliminate Φ v t from ( 18) by replacing it with the maximizer [START_REF] Ito | Crystal-MUSIC: Accurate localization of multiple sources in diffuse noise environments using crystal-shaped microphone arrays[END_REF]. Therefore, the minimization of J LS (Θ) w.r.t. Θ is equivalent to the minimization of

Φ v t = P[Φ y t ]-ϕ x t P[B x ] of
T ∑ t=1 ∥P ⊥ [Φ y t ] -ϕ x t P ⊥ [B x ]∥ 2 (19)
w.r.t. {ϕ x t } 1≤t≤T and B x . Since [START_REF] Srebro | Weighted low-rank approximations[END_REF] does not depend on P[B x ], it has infinitely many optimal solutions. This indeterminacy is resolved by exploiting the sparseness of the eigenvalues of Φ x t [START_REF] Ito | Robust Microphone Array Signal Processing against Diffuse Noise[END_REF]. In Section 4.2, we propose first estimating the V ⊥ -component, and then completing the missing V-component through matrix completion techniques. In Section 4.3, on the other hand, we propose minimizing J LS (Θ) directly subject to a rankone constraint on B x .

Optimization algorithm based on low-rank matrix completion

Based on the observation in Section 4.1, the first algorithm for least squares estimation follows the following procedure:

1. estimate Z P ⊥ [B x ] and {ϕ x t } 1≤t≤T by minimiz- ing ∑ T t=1 ∥P ⊥ [Φ y t ] -ϕ x t P ⊥ [B x ]∥ 2 , 2.
estimate B x through low-rank matrix completion [START_REF] Srebro | Weighted low-rank approximations[END_REF] of Z, 3. reestimate {ϕ x t } 1≤t≤T using B x . The algorithm eliminates {Φ v t } 1≤t≤T as a nuisance parameters, and it estimates the rest. We denote the parameter set in this method by Ω {{ϕ x t } 1≤t≤T , h}. In the first step, if we eliminate the nuisance parameter Φ v t by replacing it with its optimal value

Φ v t ← P[Φ y t ] -P[Φ x t ], (20)
the cost function reduces to

T ∑ t=1 ∥P ⊥ [Φ y t ] -ϕ x t P ⊥ [B x ]∥ 2 = T ∑ t=1 ∥P ⊥ [Φ y t ] -ϕ x t Z∥ 2 . ( 21
)
This can be minimized w.r.t. ϕ x t and Z by alternately applying the following update rules

ϕ x t ← ⟨P ⊥ [Φ y t ], Z⟩ ∥Z∥ 2 , ( 22
)
Z ← ∑ t ϕ x t P ⊥ [Φ y t ] ∑ t (ϕ x t ) 2 , ( 23
)
which is based on coordinate descent. Prior to the iterations, Z can be initialized roughly, e.g., by eigenvalue truncation of Φ y t . In the second step, we estimate B x using Z = P ⊥ [B x ] obtained in the first step, by low-rank matrix completion techniques [START_REF] Ito | Robust Microphone Array Signal Processing against Diffuse Noise[END_REF][START_REF] Srebro | Weighted low-rank approximations[END_REF][START_REF] Toh | An accelerated proximal gradient algorithm for nuclear norm regularized linear least squares problems[END_REF]. The low-rank matrix completion techniques are originally techniques for completing missing entries of a low-rank matrix, and we have extended them to the completion of a missing subspace [START_REF] Ito | Robust Microphone Array Signal Processing against Diffuse Noise[END_REF]. Here, we apply this technique to the completion of the missing matrix subspace V. As an example, we consider here the following optimization problem under a low-rank constraint:

min B x ∥P ⊥ [B x ] -Z∥ 2 (24)
s.t. B x : Hermitian positive semidefinite, rank(B x ) ≤ R.

That is, we seek for the matrix B x with a rank no larger than R (a predetermined upper bound of the rank), for which P ⊥ [B x ] is closest to Z obtained in the first step.

Other techniques for low-rank matrix completion such as trace norm minimization [START_REF] Ito | Robust Microphone Array Signal Processing against Diffuse Noise[END_REF][START_REF] Toh | An accelerated proximal gradient algorithm for nuclear norm regularized linear least squares problems[END_REF] can also be applied.

The criterion [START_REF] Yoshii | Infinite positive semidefinite tensor factorization for source separation of mixture signals[END_REF] can be optimized by alternately iterating the following update rules [START_REF] Ito | Robust Microphone Array Signal Processing against Diffuse Noise[END_REF][START_REF] Srebro | Weighted low-rank approximations[END_REF]:

1. Y ← P[B x ] + Z 2. (λ r , µ r ) ←the rth largest eigenvalue and the corre- sponding eigenvector of Y 3. B x ← ∑ R r=1 max{λ r , 0}µ r µ H r
It is guaranteed that the above procedure decreases the cost function in (24) monotonically. B x can be initialized by

B x ← ∑ t ϕ x t Φ y t ∑ t (ϕ x t ) 2 , ( 25
)
which can be computed using ϕ x t estimated in the first step. ( 25) is derived by the minimization of ∑ t ∥Φ y tϕ x t B x ∥ 2 obtained by neglecting the noise contribution in the cost function [START_REF] Zelinski | A microphone array with adaptive postfiltering for noise reduction in reverberant rooms[END_REF].

In the final step, Z = P ⊥ [B x ] is updated using B x estimated in the second step, and ϕ x t is reestimated by [START_REF] Ozerov | Multichannel nonnegative matrix factorization in convolutive mixtures for audio source separation[END_REF].

The algorithm is summarized as follows, with iter num denoting the number of iterations:

Algorithm 1

Initialize Z by rank-one approximation of Φ y t for cnt = 1 to iter num do for t = 1 to T do

ϕ x t ← ⟨P ⊥ [Φ y t ], Z⟩ ∥Z∥ 2 end for Z ← ∑ t ϕ x t P ⊥ [Φ y t ] ∑ t (ϕ x t ) 2 end for Initialize B x by B x ← ∑ t ϕ x t Φ y t ∑ t (ϕ x t ) 2 for cnt = 1 to iter num do Y ← P[B x ] + Z (λ r , µ r ) ←the rth largest eigenvalue and the corre- sponding eigenvector of Y B x ← ∑ R r=1 max{λ r , 0}µ r µ H r end for for t = 1 to T do ϕ x t ← ⟨P ⊥ [Φ y t ], P ⊥ [B x ]⟩ ∥P ⊥ [B x ]∥ 2 ϕ x t ← max{ϕ x t , 0} end for

Optimization algorithm based on a rank-one constraint

The second method minimizes the cost function [START_REF] Zelinski | A microphone array with adaptive postfiltering for noise reduction in reverberant rooms[END_REF], subject to the rank-one constraint

B x = hh H . ( 26
for cnt = 1 to iter num do for t = 1 to T do Φ v t ← P[Φ y t ] -ϕ x t P[hh H ] ϕ x t ← max { h H (Φ y t -Φ v t )h, 0 } end for h ← unit principal eigenvector of ∑ t ϕ x t (Φ y t -Φ v t ) end for for t = 1 to T do ϕ x t ← |h 1 | 2 ϕ x t end for h ← h/h 1

Spatial covariance matrix decomposition based on maximum likelihood

In this section, we propose a method for spatial covariance matrix decomposition based on maximum likelihood estimation.

= -2M log π -log det Φ x t -x H t (Φ x t ) -1 x t (28) -log det Φ v t -(y t -x t ) H (Φ v t ) -1 (y t -x t )
. Marginalizing this w.r.t. x t , we have the following marginal distribution of the observed signals y t :

p(y t ; Θ) = N C (y t ; 0, Φ x t + Φ v t ). (29)

Objective function: likelihood

We estimate the parameters by maximizing the following log-likelihood function subject to Φ v t ∈ V:

J ML (Θ) = T ∑ t=1 log p(y t ; Θ) (30) = T ∑ t=1 log N C (y t ; 0, Φ x t + Φ v t ). (31)
Here, Θ is the parameter set defined by

Θ { {ϕ x t } 1≤t≤T , B x , {Φ v t } 1≤t≤T } . ( 32
)
The maximization of the likelihood can also be viewed as the minimization of a matrix Itakura-Saito divergence [START_REF] Pham | Blind separation of instantaneous mixtures of non stationary sources[END_REF][START_REF] Ozerov | Multichannel nonnegative matrix factorization in convolutive mixtures for audio source separation[END_REF]. Indeed, J ML (Θ) can be rewritten as

J ML (Θ) = - T ∑ t=1 D IS ( Φ y t ; Φ x t + Φ v t ) + const., (33)
where D IS denotes the following matrix Itakura-Saito divergence:

D IS ( Φ y t ; Φ x t + Φ v t ) Tr { Φ y t (Φ x t + Φ v t ) -1 } (34) -log det { Φ y t (Φ x t + Φ v t ) -1 } -M. Since D IS ( kA, kB ) = D IS ( A, B
) holds for any positive scalar k, we see that D IS is suitable for audio signals, which have logarithmic characteristics. The divergence is a generalization of the Itakura-Saito divergence [START_REF] Itakura | Analysis synthesis telephony based on the maximum likelihood method[END_REF], which has proven to be effective for speech processing.

Optimization algorithm based on expectation-maximization

Since J(Θ) contains the term log det(Φ x t + Φ v t ) and (Φ x t + Φ v t ) -1 , it is difficult to obtain a close-form expression of the optimal solution. On the other hand, regarding {x t } 1≤t≤T as hidden variable, and based on the expectation-maximization (EM) algorithm [START_REF] Dempster | Maximum likelihood from incomplete data via the EM algo-rithm[END_REF], we can derive an efficient optimization algorithm. The EM algorithm iterates the following E-step and M-step alternately, and it is guaranteed to converge to a local optimum.

-E-step: Update the posterior distribution p(x t |y t ; Θ ′ ) of x t , using the current estimate of Θ ′ . -M-step: Update Θ by maximizing the Q-function

Q(Θ; Θ ′ ) = T ∑ t=1 ⟨ log p(x t , y t ; Θ) ⟩ p(xt|yt;Θ ′ ) . ( 35
)
Here, ⟨•⟩ p(xt|yt;Θ ′ ) denotes the expectation w.r.t. to the current posterior probability p(x t |y t ; Θ ′ ) of the hidden variable x t .

In the E-step, we update the posterior distribution, which can be derived as follows. From (28) and the Bayes rule,

log p(x t |y t ; Θ) c = log p(x t , y t ; Θ) (36) c
= -

( x t -µ x|y t ) H ( Φ x|y t ) -1 ( x t -µ x|y t ) . ( 37
)
Here,

µ x|y t Φ x t (Φ x t + Φ v t ) -1 y t , ( 38
)
Φ x|y t Φ x t (Φ x t + Φ v t ) -1 Φ v t , (39) and c
= means the equality up to a difference of a constant independent of x t . Therefore, the posterior probability is given by

p(x t |y t ; Θ) = N C ( x t ; µ x|y t , Φ x|y t ) . ( 40
)
Because the Gaussian distribution (40) is completely determined by its mean µ x|y t and covariance matrix Φ x|y t , in the E-step it suffices to update them.

In the M-step, we update the parameters so that the Q-function is maximized. The explicit form of Qfunction is given by the following expression: ) H end for

Q(Θ; Θ ′ ) = - T ∑ t=1 log det Φ x t (41) - T ∑ t=1 Tr [ (Φ x t ) -1 ⟨ x t x H t ⟩ ] - T ∑ t=1 log det Φ v t - T ∑ t=1 Tr [ (Φ v t ) -1 ⟨ (y t -x t )(y t -x t ) H ⟩ ] . Algorithm 3 for t = 1 to T do µ x|y t ← Φ x t (Φ x t + Φ v t ) -1 y t Φ x|y t ← Φ x t (Φ x t + Φ v t ) -1 Φ v
B x ← 1 T T ∑ t=1 1 ϕ x t Φx t for t = 1 to T do ϕ x t ← 1 M Tr [ (B x ) -1 Φx t ] Φ x t ← ϕ x t B x Φ v t ←    P [ Φv t ] , V = V uncor , V BND , V real 1 M Tr [ (B v ) -1 Φv t ] B v , V = V coh end for

Comparison between least squares and maximum likelihood estimation

Compared to the least squares estimation, the Itakura-Saito divergence in the maximum likelihood estimation has the advantage of scale invariance, and it is expected to be more suitable for audio signals, which have a logarithmic nature. Furthermore, maximum likelihood estimation is based on a generative model of the observed signal, and therefore it can be integrated with other speech enhancement techniques (e.g., source separation [START_REF] Duong | Underdetermined reverberant audio source separation using a full-rank spatial covariance model[END_REF][START_REF] Sawada | Multichannel extensions of non-negative matrix factorization with complex-valued data[END_REF][START_REF] Vincent | From blind to guided audio source separation[END_REF], dereverberation [START_REF] Nakatani | Speech dereverberation based on variance-normalized delayed linear prediction[END_REF]) based on a generative model. This is a great advantage in extending the proposed methods to a versatile method that can be applied to various real-world environments.

Experimental performance evaluation on real-world data

We experimentally validated the proposed methods on real-world data.

Experimental conditions

We fabricated a 12-channel spherical microphone array with microphones at the vertices of an icosahedron of diameter 15 cm (see Fig. 1). The microphones were mounted on a rigid spherical shell. With the array, we recorded the signal and noise images in a room at the University of Tokyo (Fig. 2). The signal image was recorded while the loudspeaker played female speech [START_REF] Kurematsu | ATR Japanese speech database as a tool of speech recognition and synthesis[END_REF], and the noise image while the windows were open. They were mixed to generate the observed signals.

We compared the following four methods for spatial covariance matrix decomposition:

conv-LS: conventional least squares method proposed in [START_REF] Ito | Diffuse noise suppression using crystal-shaped microphone arrays[END_REF], prop-LS1: proposed least squares method with lowrank matrix completion in Section 4.2, prop-LS2: proposed least squares method with a rank-1 constraint in Section 4.3, prop-ML: proposed maximum-likelihood method in Section 5.

Each method was combined with the following four noise models (see Section

Experimental results

Fig. 3 shows the spectrograms and the output SNR of the observed, enhanced, and reference signals (noise model: BND). The maximum-likelihood method ("prop-ML") has resulted in a larger SNR than the least squares methods ("prop-LS1" and "prop-LS2"), mainly because of less signal distortion. Furthermore, plotted in the logarithmic scale, the spectrogram of the ML method resembles the reference signal more, which can be explained by the logarithmic nature of the Itakura-Saito divergence. Especially, the LS methods attenuated some frequency components, while the ML method avoided this issue because the logarithmic measure strongly penalizes zeros in the estimated spectra.

Table 2 shows the output SNRs (dB) [START_REF] Ito | Diffuse noise suppression using crystal-shaped microphone arrays[END_REF] of the observed and the enhanced signals. We show results obtained with oracle initialization (calculated with the reference signal) to examine the sensitivity of the noise suppression performance to initialization. The maximumlikelihood method with the BND model yielded the highest output SNR of 8.6 dB. Comparing prop-LS1, prop-LS2, and prop-ML, prop-ML gave the highest SNR for all noise models except coh. For coh, the SNR of prop-LS2 was the highest (5.2 dB), while that of prop-ML was -0.9 dB. The algorithm prop-ML did not work well for coh. This is because the sine cardinal noise coherence matrix approaches a rank-one (hence, singular) matrix at low frequencies, and hence the update of Φ v t , which includes inversion of the matrix, diverges.

Conclusion

In this paper, we have proposed a blind method for suppressing nonstationary diffuse noise. Based on the isotropic noise models, we proposed methods for spatial covariance matrix decomposition. The decomposition algorithms are based on least squares or maximumlikelihood estimation. In the experimental evaluation, the maximum likelihood estimation has resulted in a superior performance to the least squares estimation.

Fig. 1

 1 Fig. 1 Fabricated 12-element spherical microphone array of diameter 15 cm. The microphones are mounted on a rigid spherical shell.

- 20 Fig. 2

 202 Fig.2The room layout in the experiment.

Fig. 3

 3 Fig. 3 Spectrograms of the observed, enhanced, and reference signals. The SNRs of these signals were as follows: conv-LS: 7.7 dB; prop-LS1: 4.5 dB; prop-LS2: 4.9 dB; prop-ML: 8.6 dB; observed: -0.2 dB. The plot range of the power has been equalized for all spectrograms

Table 1

 1 Detailed experimental conditions

Table 1 .

 1

	The observed signals were analyzed by the short-time
	Fourier transform (STFT). The lowest 14 frequency
	bins (below 100 Hz) were discarded, which contained
	only noise. The observed spatial covariance matrix for
	the least squares methods was computed locally by av-
	eraging y t y H t over 48 consecutive frames. The other con-
	ditions are summarized in

3):

coh: fixed noise coherence model, uncor: spatially uncorrelated noise model, -BND: blind noise decorrelation model, real: real-valued noise covariance model.

Table 2

 2 Objective evaluation of the proposed methods in terms of the SNR (dB). For the blind initialization, the highest SNR is shown in boldface, and the highest SNR within each noise model is shown in Italic. The input SNR was -1.0 dB.

	method	conv-LS		prop-LS1			prop-LS2			prop-ML
	noise model	BND	coh uncor BND real coh uncor BND real	coh	uncor	BND	real
	blind initialization	7.7	4.7	7.5	4.5	-6.5 5.2	7.5	4.9	-6.5 -0.9	8 .4	8.6	5 .8
	oracle initialization	9.2	4.8	7.5	4.5	-6.5	5.2	7.5	4.9	-6.5 -0.3	8.3	11.6	17.6

Here, we omitted a constant, and abbreviated ) ′ as follows:

The update rules for the M-step can be derived by maximizing the Q-function (41) w.r.t. the parameters (see Appendix A). The following summarizes one iteration of the algorithm derived in the above:

A Derivation of the update rules in the M-step of maximum likelihood estimation By setting the partial derivative of the Q-function (41) to zero, we have

By solving this w.r.t. ϕ x t , we get [7]

Here, we defined

Next, partial differentiation w.r.t. B x gives

Solving this w.r.t. B x , we have [7]

The update rule for Φ v t depends on the explicit form of the matrix subspace V. In the following, we first show that for the class of V satisfying

we can derive a unified update rule. Clearly, the subspaces V uncor , V BND , V real defined in Section 3 belong to the class. We then derive the update rule for V coh , which does not belong to the class. When V satisfies (52), the terms of (41) depending on Φ v t can be rewritten as

Here, P[•] denotes the orthogonal projection onto V defined using the standard inner product ⟨A, B⟩ Tr[AB] of H:

Here, {Q d } 1≤d≤D is an orthonormal basis of V, and D denotes the dimension of V. The explicit form of Q d depends on the choice of V, for which the readers are referred to [START_REF] Ito | Robust Microphone Array Signal Processing against Diffuse Noise[END_REF][START_REF] Ito | General algorithms for estimating spectrogram and transfer functions of target signal for blind suppression of diffuse noise[END_REF]. The term in P[•] in (53) generally has both components parallel and orthogonal to V. However, the latter vanishes owing to (Φ v t ) -1 ∈ V, and hence (53). To derive Φ v t ∈ V that maximizes (53), we forget the constraint Φ v t ∈ V for the moment, and differentiate (53) w.r.t. Φ v t . We have

where

As is clear from the definition of P[•], (55) certainly satisfies Φ v t ∈ V. Although we have derived (55) through partial differentiation, we can also derive it more intuitively as follows. Inverting the sign and ignoring a constant independent of Φ v t , (53) becomes the following matrix Itakura-Saito divergence:

Therefore, the maximization of (53) is equivalent to the minimization of (58 ] . Next we consider the case V = V coh . Substituting

into the Q-function (41), and differentiating it w.r.t. ϕ v t , we have, as in the derivation of (47),


Dao Trinh

Xosé Michel Keddam

Ramon Novoa

Vincent Vivier
email: vincent.vivier@upmc.fr

Xosé Ramón

Characterization of adsorbates by transient measurements in Scanning Electrochemical Microscopy

Keywords: Scanning Electrochemical Microscopy, AC-SECM, Adsorbate intermediate, Competitive adsorption, Electrochemical impedance spectroscopy

teaching and research institutions in France or abroad, or from public or private research centers.

Introduction

Electrochemical impedance spectroscopy, EIS, has been widely used for investigating complex electrochemical mechanisms such as those involved in corrosion or metal dissolution [START_REF] Keddam | Reaction model for iron dissolution studied by electrode impedance. II. Determination of the reaction model[END_REF][START_REF] Keddam | Reaction model for iron dissolution studied by electrode impedance. I. Experimental results and reaction model[END_REF][START_REF] Epelboin | Model of the anodic behavior of iron in sulfuric acid medium[END_REF][START_REF] Cachet | EIS investigation of zinc dissolution in aerated sulfate medium. Part I: bulk zinc[END_REF][START_REF] Cachet | EIS investigation of zinc dissolution in aerated sulphate medium. Part II: zinc coatings[END_REF]. When combined with one or more other potential-modulation based techniques, it allows multiple transfer functions to be determined offering the possibility of gathering simultaneously various type of data. For instance, the kinetics and the exchange of ions and molecules inside a thin polymer film has been investigated using ac-electrogravimetry method [START_REF] Gabrielli | Electrogravimetry on conducting polymers. Application to polyaniline[END_REF][START_REF] Agrisuelas | Kinetic and Mechanistic Aspects of a Poly(o-Toluidine)-Modified Gold Electrode. 2. Alternating Current Electrogravimetry Study in H2SO4 Solutions[END_REF][START_REF] To | Redox Switching of Heteropolyanions Entrapped in Polypyrrole Films Investigated by ac Electrogravimetry[END_REF], whereas the use of a double-excitation technique allowed the investigation of the double layer relaxation [START_REF] Antano-Lopez | A new experimental approach to the time-constants of electrochemical impedance: frequency response of the double layer capacitance[END_REF][START_REF] Antano-Lopez | Interface capacitance at mercury and iron electrodes in the presence of organic compound[END_REF][START_REF] Larios-Duran | Dynamics of double-layer by AC Modulation of the Interfacial Capacitance and Associated Transfer Functions[END_REF].

It has also been shown that the rotating ring-disk electrode (RRDE) allows the investigation of electrochemical and coupled chemical reactions [START_REF] Albery | Ring-disk electrodes. VII. Homogeneous and heterogeneous kinetics[END_REF][START_REF] Albery | Ring-disk electrodes. V. First-order kinetic collection efficiencies at the ring electrode[END_REF][START_REF] Albery | Ring-disk electrodes. II. Theoretical and experimental collection efficiencies[END_REF]. In their pioneering work, the group of Albery has used the alternating current measurements for defining a complex collection efficiency as the ratio of the ac-current collected at the ring to the ac-current generated at the disk [START_REF] Bruckenstein | Ring-disk electrodes. Part 17. Ring response to periodic disk electrode forcing functions[END_REF][START_REF] Albery | Ring-disk electrodes. Part 18. Collection efficiency for high frequency a.c[END_REF][START_REF] Albery | Ring-disk electrodes. Part 19. Adsorption studies at low frequency a.c[END_REF]. This approach was successfully used for characterizing the formation of thin films at an electrode surface [START_REF] Benzekri | The a.c. response of a rotating ring-disk electrode: application to 2-D and 3-D film formation in anodic processes[END_REF][START_REF] Benzekri | The a.c. response of RRDE during the passivation of iron[END_REF], and was extended to channel flow double electrode for the study of the iron dissolution [START_REF] Itagaki | Active dissolution mechanisms of iron using EIS with channel flow double electrode. Influences of chloride and fluoride ions[END_REF].

Scanning electrochemical microscopy (SECM) also allows the investigation of the electrochemical interface using a microelectrode as a collecting electrode [START_REF] Wittstock | Scanning electrochemical microscopy for direct imaging of reaction rates[END_REF]22].

SECM has been successfully used to investigate intricate systems such as those encountered in corrosion [START_REF] Jensen | Application of SECM to corrosion studies[END_REF][START_REF] Casillas | Scanning electrochemical microscopy of precursor sites for pitting corrosion on titanium[END_REF][START_REF] Casillas | Pitting Corrosion of Titanium[END_REF][START_REF] Jensen | Studies of electron transfer at aluminum alloy surfaces by scanning electrochemical microscopy[END_REF][START_REF] Gabrielli | Development of a Coupled SECM-EQCM Technique for the Study of Pitting Corrosion on Iron[END_REF][START_REF] Gabrielli | A SECM assisted EQCM study of iron pitting[END_REF][START_REF] Gabrielli | Single pit on iron generated by SECM[END_REF], for performing the dosing of adsorbed intermediates (electro)generated on a substrate using the surface interrogation mode (SI-SECM) [START_REF] Rodriguez-Lopez | Interrogation of Surfaces for the Quantification of Adsorbed Species on Electrodes: Oxygen on Gold and Platinum in Neutral Media[END_REF][START_REF] Wang | Reaction of Br2 with Adsorbed CO on Pt, Studied by the Surface Interrogation Mode of Scanning Electrochemical Microscopy[END_REF][START_REF] Rodriguez-Lopez | Reaction of various reductants with oxide films on Pt electrodes as studied by the surface interrogation mode of scanning electrochemical microscopy (SI-SECM): Possible validity of a Marcus relationship[END_REF]. Recently, Trinh et al. [START_REF] Trinh | Alternating-Current Measurements in Scanning Electrochemical Microscopy, Part 1: Principle and Theory[END_REF][START_REF] Trinh | Alternating Current Measurements in Scanning Electrochemical Microscopy, Part 2: Detection of Adsorbates[END_REF] have shown that a transient SECM (similarly to what has been done for the RRDE) allowed a quantitative study of adsorption processes.

It is generally difficult to identify intermediate species involved in multistep reactions.

For instance, halide and hydrogen adsorption have a competitive character [START_REF] Garcia-Araez | On the electrochemical behavior of the Pt(100) vicinal surfaces in bromide solutions[END_REF][START_REF] Garcia-Araez | Thermodynamic studies of chloride adsorption at the Pt(111) electrode surface from 0.1 M HClO4 Solution[END_REF][START_REF] Garcia-Araez | Competitive adsorption of hydrogen and bromide on Pt(100): Mean-field approximation vs. Monte Carlo simulations[END_REF] and various mechanisms involved in corrosion refer to the formation of adsorbates or thin films [START_REF] Keddam | Reaction model for iron dissolution studied by electrode impedance. II. Determination of the reaction model[END_REF][START_REF] Keddam | Reaction model for iron dissolution studied by electrode impedance. I. Experimental results and reaction model[END_REF][START_REF] Benzekri | The a.c. response of RRDE during the passivation of iron[END_REF].

In this paper, we report on possibility of studying elementary steps by transient SECM.

Two different examples are presented: the reduction of proton on a polycrystalline Pt electrode compared with the competitive adsorption of bromide ions and the dissolution of iron in acidic solution. In both cases, the mechanisms are assumed to involve adsorbate intermediates.

Experimental

Experimental setup

All experiments were performed with a homemade SECM device already described elsewhere [START_REF] Gabrielli | Concentration mapping around copper microelectrodes studied by scanning electrochemical microscopy[END_REF][START_REF] Gabrielli | Electrochemical Impedance Spectroscopy Investigations of a Microelectrode Behavior in a Thin-Layer Cell: Experimental and Theoretical Studies[END_REF]. Briefly, it consisted in a 3-axis positioning system (VP-25XA, Newport) driven by a motion encoder (ESP300, Newport) allowing a 3-axis displacement with a 100 nm resolution. The electrochemical measurements were performed with a homemade bi-potentiostat coupled to a low noise current-to-voltage transducer (Femto DLPCA200, BFI Optilas) with adjustable gain (10 3 -10 11 V/A). Current and potential measurements were performed with an AD/DA card (PCI-6281 -National Instruments), and the whole experimental device is controlled by a homemade software allowing usual SECM experiment to be achieved, but also multiple transfer functions to be obtained from the different input channels of the data acquisition card. A four-electrode configuration (2 working electrodes, one reference electrode, and one counter electrode) was used as described previously [START_REF] Trinh | Alternating-Current Measurements in Scanning Electrochemical Microscopy, Part 1: Principle and Theory[END_REF][START_REF] Trinh | Alternating Current Measurements in Scanning Electrochemical Microscopy, Part 2: Detection of Adsorbates[END_REF].

Materials and Reagents

Microelectrodes were made from 10 and 40 µm-diameter platinum wires, sealed in a glass capillary by using a microforge. A polycrystalline platinum electrode of 5 mm in diameter isolated by an epoxy resin was used as substrate.

All solutions were prepared in distilled water from analytical-grade chemicals. A 0.5 M H 2 SO 4 solution was used for the study of the proton reduction and the adsorption of hydrogen on polycrystalline platinum electrode. In order to study the competitive adsorption between hydrogen and bromine on platinum, a mixture of 0.1 M HClO 4 and 5 mM KBr was used. Prior to the experiments, Pt electrodes were cleaned by performing repetitive scans over one hour between the HER and the OER reaction in sulfuric acid solution.

The iron used in this study was a massive high-purity iron (99.99%) purchased from Johnson-Matthey. It was firstly annealed during 4 h at 995°C in vacuum and then quenched in water. The rod was laterally insulated with a cataphoretic paint and embedded in an epoxy resin. Prior to each experiment, the electrode was polished with SiC paper down to grade 2400, cleaned in water in an ultrasonic bath and then dried with warm air. The anodic dissolution of iron was investigated in 0.5 M H 2 SO 4 solution.

SECM measurements

The dimensionless size of the microelectrode tip RG (RG = r g /a, where r g is the total external diameter of the microelectrode and a the disc microelectrode radius) and the dimensionless tip-to-substrate distance L (L = d/a where d is the actual tip-to-subtrate distance) were determined by performing approach curves in negative feedback mode.

The SECM experiments in transient mode and the electrochemical impedance spectroscopy were performed under potentiostatic regulation of the substrate with a 50 mV peak-to-peak sine wave perturbation signal. This amplitude was chosen to improve the signal-to-noise ratio, while the linearity condition for the electrochemical system is maintained. The frequency was varied from 10 kHz down to 50 mHz with 7 frequencies per decade.

Basic of the transient SECM measurements

The transient technique used in this study is based on the former works of Albery et al. [START_REF] Albery | Ring-disk electrodes. VII. Homogeneous and heterogeneous kinetics[END_REF][START_REF] Albery | Ring-disk electrodes. V. First-order kinetic collection efficiencies at the ring electrode[END_REF][START_REF] Albery | Ring-disk electrodes. II. Theoretical and experimental collection efficiencies[END_REF] and Benzekri et al. [START_REF] Benzekri | The a.c. response of a rotating ring-disk electrode: application to 2-D and 3-D film formation in anodic processes[END_REF][START_REF] Benzekri | The a.c. response of RRDE during the passivation of iron[END_REF] devised for rotating ring-disk electrode. In generator/collector mode, an obvious analogy exists between the SECM and the RRDE:

redox species generated at one electrode can be detected at the second electrode, the timeof-flight being governed by the distance between the probe and the substrate of the SECM experiment, and by the rotation rate and the interelectrode distance in the case of the RRDE. In addition, the collection efficiency can be defined as the ratio of the current of the collecting electrode to that of the generating one. Similar to a usual electrochemical impedance measurement, a small sine-wave perturbation can be superimposed to the dc component of the substrate (at steady-state), allowing the complex transfer function   N  corresponding to the complex collection efficiency to be defined as:

      tip sub I N I       (1)
              tip tip sub tip sub sub I N I                  (2)
In the right hand-side of Eq.( 2), the first term is the inverse of the complex collection efficiency on the tip,

      tip tip tip N I      
. For this technique, a simple redox reaction has to be chosen on the tip (i.e. a kinetically fast electrochemical reaction non complicated by any coupled chemical reaction) in order to maintain its value frequency independent.

The second term is the complex transport coefficient,

      tip t sub N       , corresponding
to the mass transport contribution of electroactive species between the substrate and the tip collector electrode. This contribution always exists and can be determined experimentally. Fig. 1 shows an example of the Nyquist representation of the complex collection efficiency measured for a SECM experiment when using ferricyanide as redox mediator. The potential of the Pt substrate was set at 0.465 V/SHE and the Ptmicroelectrode (tip) was biased at 0.745 V/SHE to oxidize the collected Fe 2+ species.

When an ac-perturbation was applied to the substrate, the ferrocyanide ions diffuse in solution and are detected at the tip. The shape of the curve (a semi-cardioid) is very similar to the curves obtained with the RRDE system [START_REF] Benzekri | The a.c. response of a rotating ring-disk electrode: application to 2-D and 3-D film formation in anodic processes[END_REF][START_REF] Benzekri | The a.c. response of RRDE during the passivation of iron[END_REF] or with a channel flow double electrode system [START_REF] Itagaki | Active dissolution mechanisms of iron using EIS with channel flow double electrode. Influences of chloride and fluoride ions[END_REF]. The low frequency limit corresponds to the steady-state collection efficiency, whereas the high frequency limit was about 100 Hz, which is larger than with other generator / collector systems and depends on the probe-to-substrate distance as expected [START_REF] Itagaki | Active dissolution mechanisms of iron using EIS with channel flow double electrode. Influences of chloride and fluoride ions[END_REF]. At higher frequency, the small tail was ascribed to the ohmic coupling between the two working electrodes as already pointed out for a generator / collector system by Gabrielli et al. [START_REF] Gabrielli | Potential distribution at the surface of a ring-disk electrode[END_REF]. Thus, the collection response is similar to that of a low-pass filter [START_REF] Benzekri | The a.c. response of a rotating ring-disk electrode: application to 2-D and 3-D film formation in anodic processes[END_REF][START_REF] Itagaki | Active dissolution mechanisms of iron using EIS with channel flow double electrode. Influences of chloride and fluoride ions[END_REF] and Trinh et al. have also shown that this contribution can be evaluated numerically by using finite element simulations [START_REF] Trinh | Alternating-Current Measurements in Scanning Electrochemical Microscopy, Part 1: Principle and Theory[END_REF].

-----Figure 1 -----The last term in the right hand-side of Eq.( 2) is the complex emission efficiency of the substrate,

      sub sub sub N I       ,
sub sub sub sub d Q I n F dt      (3)
where n sub is the number of electrons involved in the electrochemical reaction at the substrate, and F the Faraday constant. In the frequency domain, Eq.( 3) can be rewritten in dimensionless form as:

    1     sub sub sub sub Q n FN j I    (4)
In a Nyquist representation Eq.( 4) corresponds to a semi circle in which

    sub sub Q j I      is
the contribution of the charge exchanged for the formation of the adsorbates at the substrate. The high frequency limit of Eq.( 4) corresponds to the kinetics of processes taking place at the substrate when controlled by the charge transfer resistance.

Conversely, at the low frequency limit, the right hand-side term is null and the steadystate requires that the flux and the current are linked through the Faraday constant.

The experimental setup used in this work allows the simultaneous measurement of the complex collection efficiency and the electrochemical impedance of the substrate,   sub Z  . Thus from Eq.( 4), the differential capacitance of the substrate, d C , can be then calculated as:

          1 sub sub sub d sub sub Q n F N C E j Z             (5)
This provides a unique way of determining the relaxation ascribed to the charge stored at the interface during the electrochemical process.

In the following section, this new technique will be used for the study of competitive adsorption between hydrogen and bromide ions on a Pt substrate and for the dissolution mechanism of iron.

Results and discussion

Competitive adsorption between Br ads and H ads

Fig. 2a shows the schematic representation of a monolayer of adsorbed hydrogen on platinum. The charge accumulated by this monolayer has been previously estimated by transient SECM at about 212 ± 12 µF cm -2 from the low frequency limit of the Nyquist plot of the differential capacitance [START_REF] Trinh | Alternating Current Measurements in Scanning Electrochemical Microscopy, Part 2: Detection of Adsorbates[END_REF].

-----Figure 2 -----In the presence of bromide ion, the schematic representation of competitive adsorption between the bromine and hydrogen is presented in Fig. 2b.

Fig. 3 shows the cyclic voltammogram of the platinum substrate in a 5 mM KBr + 0.1 M HClO 4 solution. This voltammetry curve can be used to highlight the following phenomena: the reduction of the proton at potential lower than -0.005 V/SHE, the adsorption of hydrogen in the potential range between 0.045 and 0.260 V/SHE , oxidation of the bromide at its standard potential about 1.19 V/SHE (Eq.( 6)), and adsorption of bromide ion at 0.135 V/SHE (Eq.( 7)). These values are in good agreement with literature data [START_REF] Ferro | The Bromine Electrode. Part I: Adsorption Phenomena at Polycrystalline Platinum Electrodes[END_REF][START_REF] Salaita | Structure and composition of a Platinum(111) surface as a function of pH and electrode potential in aqueous bromide solutions[END_REF] -----Figure 3 -----

2 2Br Br 2e     ( 6
) ads Br Br e     (7)
Ferro et al. [START_REF] Ferro | The Bromine Electrode. Part I: Adsorption Phenomena at Polycrystalline Platinum Electrodes[END_REF][START_REF] Ferro | The bromine electrode Part II: reaction kinetics at polycrystalline Pt[END_REF] studied the competitive adsorption between bromine and hydrogen on polycrystalline platinum electrode by gradually adding bromide ions (from 0 to 100 mM) in a solution of HClO 4 . They showed that the presence of bromide significantly alters the adsorption of hydrogen that was attributed to the competitive adsorption of bromine through the following reaction:

2 1 2 ads Br Pt H Pt Br H        (8)
The presence of bromide significantly alters the adsorption of hydrogen resulting in the gradual disappearance of the adsorption/desorption peak at about 0.260 V/SHE.

Concomitantly, the peak at 0.135 V/SHE increases with the bromide ion concentration, which confirms the competitive adsorption of bromine and hydrogen.

In this case, the total exchanged charge involved through the competitive formation of the hydrogen and bromine absorbed layer is the sum of the charge required for the adsorption of each species,

H sub Q  and Br sub Q  , respectively. H B r sub sub sub Q Q Q    (9)
The potential corresponding to the adsorption reaction of bromine is very close to the adsorption potential of hydrogen and it is difficult to distinguish between these two phenomena by conventional electrochemical methods, especially on a single voltammogram. By performing the transient SECM measurements while detecting different species generated from the substrate, specific information about each adsorption process can be gathered and the competitive reaction between the adsorption of bromine and hydrogen can be selectively studied. -----Figure 4 -----At 1.19 V/SHE, Br -is oxidized to Br 2 (Br 3 -) which is then detected at the tip (tip potential is maintained at 0.645 V/SHE to reduce the Br 2 generated at the substrate). In both cases, the HF limit tends towards 0 because in the frequency domain, the impedance is governed purely by the charge transfer resistance. The low frequency limit corresponds to the charge accumulation at the electrode surface, i.e., 1 mFcm -2 for hydrogen and 1.2 mF mFcm -2 for bromine (Table 1). The latter value corresponds to the accumulation of charge on the substrate associated with the competitive adsorption between ads H and ads Br .

-----Table 1 -----

Thus, varying the electrode potential allows the different contribution of absorbed species to be evaluated and the sum of charges involved in the formation of ads H and ads Br is

given by

H B r sub sub sub Q Q Q   
= 71 + 85 = 156 µC cm -2 . This value (156 µC cm -2 ) is smaller than the one obtained for a monolayer of adsorbed hydrogen in the absence of bromide ions (212 µC cm -2 ) which confirms the hypothesis of competitive adsorption of bromide ions on the site occupied by hydrogen on platinum substrate. It should also be mentioned that because both electrodes (tip and substrate) are made from the same material and were prepared following the same experimental procedure (grinding, washing, electrochemical cleaning), the roughness factor for each electrode was assumed to be the same. This allows expecting that the complex collection efficiency to be independent of this roughness factor.

In addition, the time-constant for the relaxation of ads H at the electrode surface without the bromide ion in solution is about 2 Hz [START_REF] Trinh | Alternating Current Measurements in Scanning Electrochemical Microscopy, Part 2: Detection of Adsorbates[END_REF]. This frequency points to a kinetics ten times faster than the value determined in presence of bromide ions (Figure 4). All the characteristic constants extracted from these experiments are reported in Table 1. This behavior indicates a change in the mechanism and kinetics of adsorption of hydrogen in presence of a competitive adsorption.

Mechanism of iron dissolution in acidic solution

The dissolution mechanism of iron in an acidic solution has been the topic of controversy [START_REF] Keddam | Mechanism of anodic dissolution of iron in acid solution[END_REF]. Two different mechanisms were considered based on the interpretation of current potential curves and electrochemical impedance diagrams [START_REF] Keddam | Mechanism of anodic dissolution of iron in acid solution[END_REF][START_REF] Bockris | The electrode kinetics of the deposition and dissolution of iron[END_REF][START_REF] Heusler | Der Einfluß der Wasserstof fionenkonzentration auf das elektrochemische Verhalten des aktiven Eisens in sauren Lösungen[END_REF], in which adsorbed monovalent iron species act as an intermediate in the transfer of the first and second electrons of Fe (mechanism proposed by J. O'M. Bockris, et al. [START_REF] Bockris | The electrode kinetics of the deposition and dissolution of iron[END_REF]) whereas Heusler [START_REF] Heusler | Der Einfluß der Wasserstof fionenkonzentration auf das elektrochemische Verhalten des aktiven Eisens in sauren Lösungen[END_REF] suggested that it acts as a catalyst in the dissolution of divalent Fe cations in acidic solutions.

These two mechanisms can be summarized as follows:

Bockris' mechanism [START_REF] Bockris | The electrode kinetics of the deposition and dissolution of iron[END_REF] is described by two successive and irreversible elementary steps involving

I ads

Fe as adsorbate intermediate:

1 k I ads Fe Fe e     ( 10
) 2 k I I I ads sol Fe Fe e     (11)
Thus, mass and charge balances can be expressed as

  1 2 1 d k k dt        (12)     1 2 1 Sub i F k k      (13)
where  is the fractional surface coverage of

I ads

Fe and  is the maximum of surface concentration of adsorbate intermediate. The flux of electrogenerated species at the substrate is then given by

2 Sub k    (14)
At steady-state, the fractional surface coverage is given by

1 1 2 k k k    (15)
The impedance of the substrate is given by

      1 2 2 1 2 2 1 2 1 2 1 sub b b k F b b k k k Z k k j                  (16)
allowing the collection efficiency to be calculated as

  1 2 2 2 2 2 1 2 sub sub b b k N bk k Z k k j                 (17)
The mechanism proposed by Heusler [START_REF] Heusler | Der Einfluß der Wasserstof fionenkonzentration auf das elektrochemische Verhalten des aktiven Eisens in sauren Lösungen[END_REF]

The mass and charge balances can then be expressed as

  1 1 1 d k k dt         (20)     1 1 2 1 2 Sub i F k k k         (21)
The flux of electrogenerated species at the substrate is also given by

2 sub k    (22)
At steady-state, the fractional surface coverage is given by

1 1 1 k k k     (23)
The impedance of the substrate is given by

      1 1 1 1 1 1 2 2 1 1 2 1 1 1 2 2 sub k b b F b b k k b k k k Z k k j                           (24)
allowing the collection efficiency to be calculated as

  1 1 1 2 2 2 1 1 sub sub k b b N bk k Z k k j                    (25)
-----Figure 5 ---- Assuming that electrochemical steps follow Tafel kinetics and using the linearized expression of the currents given by Eq.( 13) and Eq.( 21), it is possible to calculate the impedance response as shown in Fig. 6a for the Bockris' mechanism and Fig. 6b Heusler's mechanism. Independently of the kinetic constants used for the simulation, in both cases the shape of the impedance diagrams are very similar to those obtained experimentally. Both the amplitude and time-constants vary in the same way. It is thus difficult to decide which mechanism is the most suited.

-----Figure 6 ----------Figure 7 -----Fig. 7 shows the calculated collection efficiency of the substrate using the same set of parameters as those used for the impedance diagrams presented in Fig. 6. Interestingly, a different behavior is observed. In the case of Bockris' mechanism, the collection efficiency locus overlap in amplitude and a small frequency shift can be observed (Fig. 7a) whereas for Heusler's mechanism, the amplitude of the complex collection efficiency is potential dependent (Fig. 7b). As the potential increases, the high frequency limit increases corresponding to a change of the kinetic process at the substrate under the assumption of vanishing variations of the superficial concentration.

-----Figure 8 -----Fig. 8 shows the experimental complex collection efficiency of the substrate measured simultaneously with the impedance diagrams presented in Fig. 5. The shape of the transfer function does not depend on the potential and a small frequency shift is observed.

These results compare favorably with the theoretical curves presented in Fig. 7a indicating that the Bockris' mechanism allows a quantitative description of the iron dissolution valid at low current densities and low pH [START_REF] Keddam | Reaction model for iron dissolution studied by electrode impedance. II. Determination of the reaction model[END_REF][START_REF] Keddam | Reaction model for iron dissolution studied by electrode impedance. I. Experimental results and reaction model[END_REF]. Thus the transient SECM allows distinguishing between different possible mechanisms. It should be mentioned that for experiments performed at higher potentials (thus for larger steady-state dissolution current), no reliable results could have been obtained. This was ascribed to the fact that under high dissolution rates the distance between the Pt-probe and the iron substrate varies too quickly compared to the time-scale of the impedance acquisition rate.

Conclusions

In this work, we have highlighted different intakes that can be obtained using the SECM in transient mode in order to characterize an electrified interface.

In the case of competitive adsorption between halide ions and hydrogen, it was shown that it is possible to discriminate between the different species and to evaluate independently the charge exchanged for each adsorbate.

In the case of the formation of an adsorbate intermediates like Fe(I) species, the complex collection efficiency of the substrate can decide between different mechanisms and appears to be a complementary approach to the usual impedance measurements.

The perspectives of this work are the investigation of reaction mechanism on well oriented materials such as the absorption and the competitive adsorption on single crystals and the investigation of the passive film formation and passive film reactivity as a function of the grain orientation.

Figure captions

  -components of the tip and substrate currents, respectively. The introduction of the electrochemical fluxes at the substrate, .(1) allows the rewriting of   N  as:

Figure 4

 4 Figure4shows the differential capacitance

-Fig. 5

 5 Fig. 5 shows the experimental impedance diagrams measured at the onset of the dissolution peak of pure iron as a function of the electrode potential. In each case two loops are evidenced, a high frequency capacitive one corresponding to the double layer capacity-charge transfer resistance and a low frequency inductive one generated by the frequency response of the surface coverage by the adsorbed reaction intermediate. These diagrams are characterized by a decrease of the impedance (larger value of the dissolution current) and a shift of the characteristic frequency towards shorter time-constant as the potential increases for both loops. Assuming that electrochemical steps follow Tafel

Figure 1 : 6 Fe(CN) + 0. 5 MFigure 2 :

 1652 Figure 1: Complex experimental collection efficiency for a 50 mM

Figure 3 :

 3 Figure 3: Cyclic voltammogram of the polycrystalline platinum electrode (0.5 cm in diameter) in a 5 mM KBr + 0.1 M HClO 4 solution -initial potential 0.445 V/SHE -5 cycles at 50 mVs -1 .

Figure 4 :

 4 Figure 4: Differential capacitances calculated from the complex experimental collection efficiency according to Eq. (5) sub E = -0.005 V/SHE, Pt-tip: 20 μm in radius, Pt-substrate

Figure 5 :

 5 Figure 5: Electrochemical impedance diagrams of pure iron in sulfuric acid solution (0.5 M) as a function of the electrode potential.

Figure 6 :

 6 Figure 6: Simulated electrochemical impedance diagrams for pure iron in sulfuric acid solution as a function of the electrode potential using Refs. [1, 2] for the choice of kinetic

Figure 7 :

 7 Figure 7: Simulated collection efficiency of the substrate for pure iron in sulfuric acid solution as a function of the electrode potential (same kinetic constants as in Fig. 6). (a) Bockris' mechanism, (b) Heusler's mechanism.

Figure 8 :

 8 Figure 8: Experimental collection efficiency of the substrate for pure iron in sulfuric acid solution (0.5 M) as a function of the electrode potential.

Table 1 : 1 Figure 1 - 5

 1115 Figure 1

Figure 8 9

 8 Figure 8

  also involves

									I ads Fe as adsorbate intermediate,
	but the later acts as a catalyst:
	Fe	1 1 k k   	I ads Fe		e			(18)
	I ads Fe Fe 	2   k	I ads Fe		I I sol Fe		2	e 

Table 1

 1

		Capacitance	Charge	Characteristic
		/ mF cm -2	/ µC cm -2	Frequency / Hz
	Only H ads (without the presence of Br)	3	212	2
	H ads competitive with Br ads	1	71	0.2
	Br ads competitive with H ads	1.2	85	0.5

Acknowledgements

The authors gratefully acknowledge D. Rose (UPR15) for technical support.


A Nassiopoulos
email: alexandre.nassiopoulos@ifsttar.fr

J Brouns

N Artiges

M Smail

B Azerou

RETROFIT: A SOFTWARE TO SOLVE OPTIMIZATION AND IDENTIFICATION PROBLEMS APPLIED TO BUILDING ENERGY MANAGEMENT

Keywords: model calibration, state-parameter identification, inverse problems

Problems such as parameter identification for model calibration, optimal design or optimal energy management can all be formulated in a similar framework as problems consisting in finding the minimum of a cost function. The paper presents the software ReTrofiT that specifically treats this kind of problems applied to building energy performance models. ReTrofiT is first of all a simulation tool for evaluating building thermal behaviour and computing energy consumptions. The novelty compared to state-of-the-art energy simulation software is that it also integrates a generic set of tools and algorithms to set up and solve optimization problems related to the building thermal model. The use of the adjoint model, that is intrinsically implemented in the code, constructs fast and efficient algorithms to solve linear, non linear, constrained or unconstrained problems addressing a wide range of applications.

INTRODUCTION

Energy management systems in buildings greatly contribute to the improvement of overall energy efficiency. Monitoring systems can lead to significant reductions of the global energy use by increasing occupants' awareness of the consumptions or by enabling the implementation of more efficient regulation strategies. Model predictive control consists in computing optimal heating or cooling strategies by taking into account the future evolution of the state of the building under forecast weather or use conditions. Demand response strategies in smart grids consist in adjusting energy demand at the enduser level to reduce the overall demand thus resulting in end-user customer bill savings, increase of electricity market stability and of electricity supply reliability. Further, today the building construction practices evolve towards a more performance-based approach in which the concern becomes the performance of the final building rather than the means employed to construct it.

All the aforementionned applications rely on the ability to accurately predict a building's behavior using a calibrated model. In building energy applications, uncertainties in input data of modelling tools often lead to important discrepancies between the model predictions and the real performance. The desired model response can be obtained if the internal parameters of the model are calibrated using on-site measurements and model identification methods.

This kind of problems can typically be formulated as optimization problems consisting in minimizing a cost function measuring the discrepancy between sensor data and model's response [START_REF] Reddy | Literature review on calibration of building energy simulation programs: Uses, problems, procedures, uncertainty and tools[END_REF] [2], [START_REF] Heo | Calibration of building energy models for retrofit analysis under uncertainty[END_REF], [START_REF] Bacher | Identifying suitable models for the heat dynamics of buildings[END_REF]. Such optimization problems also come up in energy management applications, model reduction techniques or optimal design.

The paper presents the software ReTrofiT that was specifically designed to treat this kind of optimization problems in the context of building energy simulation. ReTrofiT is first of all a dynamic building simulation code that compares to state-of-the-art simulation software such as Energy Plus or TRNSYS. ReTrofiT adopts the so-called multizone modelling assumption for the description of the building's envelope that leads to solve a set of ordinary and partial differential equations.

The major force of ReTrofiT comes from the fact that it integrates a set of tools and algorithms to set up and solve optimization problems related to the building energy model. Together with the model used for simulation it implements the adjoint model. The adjoint model computes in an efficient way the derivatives with respect to any set of scalar or functional variables and to construct fast optimization algorithms [START_REF] Özisik | Inverse Heat Transfer[END_REF]. It thus enables to construct fast minimization algorithms and to compute sensitivities with respect to any of the model's variables. ReTrofiT also integrates regularization tools in order to deal with the ill-posed nature of some identification problems [START_REF] Engl | Regularization of Inverse Problems[END_REF]. The set of algorithms provided can tackle both linear and non-linear optimization problems, with or without constraints.

ReTrofiT is implemented as a Matlab R toolbox [START_REF]Matlab software[END_REF]. Its object-oriented architecture provides the necessary modularity so that new model of algorithm components be easily added. Any CAD tool compatible with the standard gbXML data exchange format can be used to generate the building's geometrical description using the specific import tool.

This paper presents the main characteristics of ReTrofiT software. The first section describes the building energy model. The next session describes the tools provided to set up and solve optimization problems. Finally, the last session illustrates the use of the software in typical applications related to heating load computations, sensitivity analysis and identification of non linear model parameters.

MULTIZONE MODEL AND DYNAMIC THERMAL SIMULATION

Modelling assumptions

The mathematical model describing the thermal state of the building builds upon standard multizone assumptions (see for instance [START_REF] Clarke | Energy simulation in building design[END_REF] for details). The zone state variables are supposed homogeneous in the entire air volume. In each zone, the mean air temperature is governed by an ordinary differential equation with source terms related to internal gains, solar gains, and air mass exchange gains.

Air flows are considered in a simplified way, only air mass exchange between zones and with outside air due to infiltration effects are described. Expression Definition R jk (T k -T j ) Air mass exchange between zone j and zone k; R jk is an exchange coefficient, the flow direction being set by the sign of the temperature difference (T k -T j ) R e j (T e -T j ) Air mass exchange between zone j and the outside Γ d j Solar gain coefficient (beam component) for zone j Γ f j Solar gain coefficient (diffuse component) for zone j Short wave solar radiation gains are taken into account through Γ d j and Γ f j . They are computed from the number of windows in the zone, the relative orientation and the window transmission coefficient.

       C j dT j dt -∑ i C 0 i j -∑ i C L i j = R jk (T k -T j ) + R e j (T e -T j ) +Γ d j Φ d + Γ f j Φ f + Q j +W j + η j P T j (0) = 0 (1)
                           c i ∂ θ i ∂t - ∂ ∂ x k i ∂ θ i ∂ x = f i [0, L i ] × [0, τ] -k i ∂ θ i (0,t) ∂ x = ∑ j C 0 i j + C 0 i0 + ∑ k,m R 0m ik + γ d i0 Φ d + γ f i0 Φ f [0, τ] k i ∂ θ i (L i ,t) ∂ x = ∑ j C L i i j + C L i i0 + ∑ k,m R L i m ik + γ d iL i Φ d + γ f iL i Φ f [0, τ] θ i (x, 0) = 0 [0, L i ] (2)
In the above expressions, x n and x m takes the values 0 or L i depending of the surfaceorientation. The table below gives the definitions of the various terms.

Expression

Definition

C n i j = h 0 i j (θ i (x n ) -T j )
Convection terms between wall i (at x = x n ) and zone j at temperture

T j C n i0 = h n i0 (θ i (x n ) -T e )
Convection terms between wall i (at x = x n ) and outside air at temperature

T e R nm ik = α nm ik (θ i (x n ) -θ k (x m )) Radiative transfer between wall i (at x = x n ) and wall k (at x = x m ) R n∞ i j = β 0 i j (θ i (x n ) -T ∞ )
Radiative transfer between wall i (at x = x n ) and the outside environment at temperature T ∞ γ d in , γ f in Solar coefficient of face x = x n of wall i, for the beam et diffuse radiation

Heating device dynamics are governed by ordinary differential equation describing the evolution of the amount of heat W j supplied into zone j when some energy power P j is injected.

   d j dW j dt +W j = η j P j W j (0) = W 0 j (3)
This generic model can represent, in a simplified way, both systems composed of standalone devices (eletric convectors for instance) and centralized heat production facilities. Equations ( 1) to (3) constitute the standard system of equations in a typical energy simulation. As mentioned above, all equation terms are implemented as independent features that can be activated or disabled. These base equations can lead to different levels of details in the building model. Further, the software architecture enables future extensions of the system of equations with addition of new terms and equations, probably from other engineering fields.

Software architecture

ReTrofiT software is built in a modular way with object oriented programming. The building model constituting the core of the software is defined by a class containing all functions for mesh generation, matrix assembly and numerical integration. Each equation of the model is described into a separate object, which contains the implementation of all routines needed for the computation of its own terms. Those objects, called features, can be disabled on will, making it easy to adapt the level of detail.

In addition, this particular architecture makes it easy to add new terms in existing equations or new equations in the whole system of equations. The software can thus address problems related to other engineering fields and could be viewed as a generic optimization platform for multiphysics applications, as the next section describes.

SOLVING OPTIMIZATION PROBLEMS

General formulation

It is a standard approach to formulate model calibration as a parameter identification problem [START_REF] Brouns | Identification de paramètres et séparation de sources thermiques à dynamiques différentes dans le bâtiment en utilisant la variation totale[END_REF]. Usually, a preliminary sensitivity analysis determines the model variables that have the highest impact on the computation of the quantities of interest. The parameter identification problem is an inverse problem that aims to determine the values of these predominant variables so that the gap between model simulation results and measured data be minimized. Such a problem is often ill-posed and needs an appropriate regularization procedure to be solved.

Problems of similar nature are also encountered in continuous monitoring applications: in this case, the unknown variables can be the initial state and/or the solicitations received by the system. The state estimation problem consists in identifying these unknowns using sensor data in order to reconstruct the global state of the system.

Optimal control for energy management applications provides another example of an optimization problem. In its simplest way, an optimal control problem consists in computing control laws (for example, power set-point of a heating device) to drive one or several state variables to a given set-point [START_REF] Hazyuk | Optimal temperature control of intermittently heated buildings using model predictive control: Part ii -control algorithm[END_REF], [START_REF] Oldewurtel | Use of model predictive control and weather forecasts for energy efficient building climate control[END_REF].

All these problems can be formulated as minimization problems : one looks for a vector u such that min

u∈U c J(u) (4)
where u is the control vector, U c the control space that takes in account any additional constraints on u. J is the cost functional which in a generic form writes

J(u) = α 1 2 Tu -y 2 M + α 2 2 u -u 0 2 U + α 3 2 φ (u) 2 (5)
The first term of this expression is the so-called residue term. It measures the quadratic gap between the model response Tu and sensor data (or set-points, depending on the problem type) in the measurement space M . The model response Tu is obtained by applying an observation operator that models the measurement process on the solution of the partial derivatives equation system (1) - [START_REF] Heo | Calibration of building energy models for retrofit analysis under uncertainty[END_REF].

The second term is a weighting term for the control vector in a specific norm. In an identification problem, it acts like a regularization term avoiding the solution to be too far from an a priori guess . In this case, the weighting coefficient α 2 is called the Tikhonov parameter and must be small [START_REF] Engl | Regularization of Inverse Problems[END_REF]. In an optimal control problem, this term represents a measure of the energy cost. In this case, coefficients α 1 and α 2 reflects the trade-off between, for instance, the comfort represented by the first term and the consumption represented by the second one.

The third term is a so-called penalization term enable to force the constraints. The function φ (u) is constructed such that it is (almost) null when u satisfies the constraints and positive otherwise.

Adjoint model

The minimum of the cost function J corresponds to the Euler conditions, i.e. to the point where the gradient is null. This minimum can be obtained using descent iterative methods.

The computation of the gradient of the second term is straightforward. The gradient corresponding to the computation of the third term depends on the function φ . In most cases, it is possible to choose a function for which the gradient is known explicitly.

Computing the gradient calculation of the first term is less straightforward because the unknowns u does not appear explicitly in the expression. Numerical differentiation should be avoided, because discretization of time dependent components of vector of unknowns u usually leads to a discrete vector of unknowns with a great number of components. Numerical differentiation would require at least one model resolution per discrete vector component which would dramatically increase computation cost.

ReTrofiT natively implements the so-called adjoint model, which gives the exact gradient with only one adjoint resolution according to optimal control theory [START_REF] Lions | Contrôle optimal des systèmes gouvernés par des équations aux dérivées partielles[END_REF]. The adjoint model has the same structure as the direct model, meaning that the same numerical tools can be used for its resolution.

Minimization algorithms

Iterative descent algorithms of descents can be used to solve (4), each iteration requiring a resolution of the adjoint models to obtain the new gradient computation. ReTrofiT implements various descent algorithms exemplified in what follows.

When the model is linear with respect to all unknowns one can use the steepest descent or conjugate gradient algorithm. The latter, applied to a functional of the form J(u) = 1 2 Tuy 2 M writes as follows :

• Initialization : u 0 • g 0 = T * (Tu 0 -y), h 0 = -g 0 • At each iteration n ≥ 1 -Compute optimal step : ρ = (g n ,h n ) h n ,T * T h n ) -Compute descent direction h n : u n+1 = u n + ρh n -Iteratively compute the gradient : g n+1 = g n + ρT * T h n -Compute conjugate direction coefficient : γ = (g n+1 ,g n+1 ) (g n ,g n )
-Iteratively compute the descent direction :

h n+1 = -g n+1 + γh n -Stop if (g n+1 , g n+1 ) < ε
In case the model is non-linear with respect to some unknowns, a linearization step is required. The Levenberg-Marquardt algorithm is suited to this case. It consists in constructing, at each iteration, the following functional:

J k (δ u) = 1 2 Tu k + δ T (δ u) -y 2 M + ε n 2 δ u U (6)
and to minimize with respect to δ u. In the expression (6), δ T (δ u) is the response of the linearized model around the point u k such that: T (u k + δ u) ∼ T (u k ) + δ T (δ u). The minimum of J k can for instance be obtained with the conjugate gradient described above. These algorithms can be adapted in the cases where bound constraints must be prescribed on the unknowns. In this case, an additional projection step is added at each iteration. For more general constraints more appropriate algorithms such as the Uzawa method or other dual algorithms should be used.

ReTrofiT software includes specific classes for all the above mentioned algorithms. Here again, the architecture is open to future addition of new algorithms.

APPLICATION EXAMPLES

This section presents application examples on a test case in order to illustrate the various possibilities of the software. The case study is a building consisting of two zones, one of which includes an opening.

Computing heating loads and sensitivities

The theoretical heating loads during a time period [0, τ] can be computed by determining W that minimizes the cost function J(W ) given by

J(W ) = 1 2 ∑ j τ 0 (T j -T c ) 2 dt + ε 2 ∑ j τ 0 W 2 j dt
where T j are the solutions of (1). The sensitivity S W of W * with respect to a parameter (for example, a perturbation δ R e j of the ir exchange rates R e j ) is given by Find S W solution of min S W J(S W ), where

J (S W ) = 1 2 ∑ j τ 0 ( Tj ) 2 dt + ε 2 ∑ j τ 0 (W * j ) 2 dt        C j dT j dt -∑ i C 0 i j -∑ i C L i j = R jk (T k -T j ) + δ R e j (T e -T j )
T (0) = 0 [START_REF]Matlab software[END_REF] where the last equation is coupled to the linearized forms of equations ( 2) and (3). The figure [START_REF] Nassiopoulos | Calibration of building thermal models using an optimal control approach[END_REF] gives an example of solutions of these two problems, obtained with ReTrofiT.

Resolution of an identification problem

The same tools can be used to formulation and solve an identification problem. Figure [START_REF] Heo | Calibration of building energy models for retrofit analysis under uncertainty[END_REF] gives an example where ambient air temperature measurements are used to reconstruct time-dependent air mass exchange. This result was obtained using the Levenberg-Marquardt algorithm with bound constraints were prescribed on the variable. The total computational time is around 11 seconds on a standard PC.

CONCLUSION

ReTrofiT is a software for building energy simulation specially designed to solve optimization and identification problems applied to building energy management. It natively implements the adjoint model and includes various algorithm classes thus providing the user with a wide set of operational tools to formulate and solve optimization problems covering a wide range of engineering applications.

Figure 1 :

 1 Figure 1 : The case study building geometry: 3D geometry description using GSketchup CAD software

Figure 2 :

 2 Figure 2 : Computation of heating load and its sensivity with respect to a perturbation on air filtration rates.

  (a) Ambient air temperature measurements (b) Reconstructed of air mass exchange

Figure 3 :Figure 4 :

 34 Figure 3 : Example of reconstruction of air mass exchange based on ambient air temperature measurements

Table 1 :

 1 Description of terms

Table 2 :

 2 Description of terms

EWSHM 2014 -Nantes, France

The software architecture was designed to adapt to various situations by enabling the user to adapt the model's level of detail and model's complexity at will. It also enables future extensions of the software possibilities by addition of new resolution algorithms or model features, possibly related to other fields of physics.