Albert Bijaoui
email: albert.bijaoui@oca.eu

A Greedy Algorithm for a Sparse Scalet Decomposition

Sparse decompositions were mainly developed to optimize the signal or the image compression. The sparsity was first obtained by a coefficient thresholding. The matching pursuit (MP) algorithms were implemented to extract the optimal patterns from a given dictionary. They carried out a new insight on the sparse representations. In this communication, this way is followed. It takes into account the goal to obtain a sparse multiscale decomposition with the different constraints: i/ to get a sparse representation with patterns looking like to Gaussian functions, ii/ to be able to decompose into patterns with only positive amplitudes, iii/ to get a representation from a translated and dilated pattern, iv/ to constrain the representation by a threshold, v/ to separate the sparse signal from a smooth baseline. Different greedy algorithms were built from the use of redundant wavelet transforms (pyramidal and à trous ones), for 1D signals and 2D images.

Experimentations on astronomical images allow one a gain of about two in sparsity compared to a classical DWT thresholding. A fine denoising is obtained. The results do not display any wavy artifacts. This decomposition is an efficient tool for astronomical image analysis.

I. INTRODUCTION

Many kinds of sparse decompositions were developed these two last decades mainly in order to optimize the signal or the image compression [START_REF] Donoho | Compressed sensing[END_REF]. Among them the wavelet transform appeared as a very efficient way to carry out sparse representations for a large class of signals [START_REF] Usevitch | A tutorial on modern lossy wavelet image compression: foundations of JPEG[END_REF]. Since different new transforms were proposed for optimizing sparsity for a given signal class [START_REF]Beyond wavelets[END_REF]. The sparsity was first obtained with an available coefficient thresholding. The Discrete Wavelet Transform (DWT) concentrates the significant information onto the minimum of coefficients without any redundancy between the data and the transform.

Sparsity can be also derived from overcomplete representations, by choosing optimally the projection vectors [START_REF] Donoho | Stable recovery of sparse overcomplete representations in the presence of noise[END_REF]. Matching pursuit (MP) algorithms were implemented to extract the representation elements from a given dictionary [START_REF] Mallat | Matching pursuits with time-frequency dictionaries[END_REF]. The nonorthogonality between the alphabet patterns can lead to convergence problems. The orthogonal matching pursuit (OMP) algorithm [START_REF] Pati | Orthogonal Matching Pursuit: Recursive Function Approximation with Applications to Wavelet Decomposition[END_REF] allows one to keep the orthogonality between the representation and the residuals. Different algorithms were proposed in order to extract in parallel the patterns [START_REF] Feichtinger | Hierarchical parallel matching pursuit[END_REF], [START_REF] Donoho | Sparse Solution of Underdetermined Linear Equations by Stagewise Orthogonal Matching Pursuit[END_REF].

The basis pursuit [START_REF] Chen | Atomic Decomposition by Basis Pursuit[END_REF] has given a new insight on the sparse representations. The introduction of a ℓ 1 constraint on the representation has reduced significantly the number of coefficients. Some algorithms were proposed in order to solve the associated variational problem. The relation between the different ways to obtain a sparse decomposition has been carefully studied. If the sparsity is sufficient, the basis pursuit algorithms minimizes also the norm ℓ 0 [START_REF] Tropp | Greed is good: algorithmic results for sparse approximation[END_REF] [START_REF] Gribonval | Beyond sparsity: Recovering structured representations by l1l1 minimization and greedy algorithms[END_REF] Here another insight is given, taking into account our goal, obtaining a sparse multiscale decomposition of large images with the following constraints:

• The dictionary has to be built from well-defined patterns.

Our goal does not consist in the search of the largest sparsity but in the search of how to get a representation with simple patterns. This is required for a further signal/image analysis for which a decomposition into separate objects is needed. The alphabet is made of monomodal, positive and symmetric patterns, looking like to Gaussian functions.

• The patterns were chosen compact. This allows one to increase the separation between them and thus improve the pattern extraction.

• The pattern amplitudes can all be positive. This constraint is connected to our scientific goal, i.e. to describe a signal/image into positive separate objects. In specific cases these amplitudes could all be negative.

• Our goal is also to get a representation as covariant as possible with translations and dilations. This leads one to generate the pattern set by translations and dilations of a basic pattern.

• The decomposition is constrained by a threshold. In the case of a noisy signal, this threshold is linked to false alarm detections at each scale.

• The signal representation is the sum of the sparse decomposition with a baseline. This baseline is such that no significant pattern can be detected from it at each examined scale, for the given threshold.

This last constraint plays an essential role in the algorithm. This characteristic is introduced in order to take into account our goal, the extraction of objects in astronomical images. In this framework, a sky background is supposed to be superimposed on the objects field. The background is first computed and then removed. Here, the background is considered as a spurious component which is simultaneously determined with the searched components.

II. THE GREEDY ALGORITHMS.

Hereafter, the algorithm is presented only in its 1D pyramidal version.

A. The basic patterns.

The multiresolution pyramid is defined by [START_REF] Mallat | A Theory for Multiresolution Signal Decomposition: The Wavelet Representation[END_REF]:

F (i, k) = 1 2 i < f (x), φ( x 2 i -k) > . (1)
where the scaling function φ(x) satisfies the dilation relation [START_REF] Strang | Wavelets and Dilation Equations: A Brief Introduction[END_REF]:

1 2 φ( x 2 ) = n=n1,n2 h(n)φ(x -n), (2)
Here h(n) is symmetric and

n 1 = -n 2 .
The following recursive relation is easily derived:

F (i + 1, k) = n=n1,n2 h(n)F (i, 2k + n). (3)
The coefficients can be also written as:

F (i, k) = m=m1,m2 h(i, m)F (2 i k + m). (4)
The patterns h(i, m) are obtained by the following recursive relation:

h(i + 1, m) = n=n1,n2 h(i, m -2 i n)h(n). (5)
At scale 0 we have h(0, m) = δ(m), which leads to h(1, n) = h(n). The bounds (m 1 , m 2 ) increase exponentially (m 1 = n 1 2 i and m 2 = n 2 2 i ). The pyramidal algorithm carries out the set of correlations of the signal with the filters {h(i, m)}. These patterns are built by a pyramidal rule (Equation 5), so, we call them pyrels for (pyramid elements). These pyrels may also be called scalet pyrels taking into account their construction from the scaling function. It can be denoted that the patterns {h(i, m)} derives from the scaling function, but it does not correspond to its sampled values.

B. The background removal.

Let us consider the wavelet transform which results from the differences between two successive approximations:

w(i + 1, k) = F (i, k) -F (i + 1, k) (6)
where F (i + 1, k) is the approximation at scale i + 1 before the decimation. The wavelet coefficients are independent of the constant adding. Taking into account the relation (4) we can write as:

w(i + 1, k) = m=m1,m2 g(i + 1, m)F (2 i k + m) (7)
where:

g(i + 1, m) = m=m1,m2 h(i, m) -h(i + 1, m) (8)
Thus, the correlations in each point between the signal and the patterns g(i, m) are obtained with this wavelet transform.

The highest amplitude coefficients can be identified and they can be used for the representation in the framework of a matching pursuit algorithm. This is equivalent to threshold the wavelet transform and to restore the signal by inversion. The result displays generally wavy artifacts. An iterative algorithm can be applied on the residuals to reduce them, until we get non significant wavelet coefficients. This process reduces the sparsity and it carries out negative artifacts without physical meaning. In order to avoid this difficulty, the signal will be restored using pyrel patterns, even if the identification is done from the wavelet coefficients. The wavelet transform will only play the role to identify the significant pyrels in presence of a background.

C. The pyrel identification.

The pyrels are selected on a criterion based on the highest SNR. In the hypothesis of a white Gaussian noise, Relation [START_REF] Feichtinger | Hierarchical parallel matching pursuit[END_REF] allows one to estimate the w(i, k) standard deviation due to the noise and thus its SNR. The wavelet coefficient w(i 0 , k 0 ) having the highest SNR identifies a pyrel at (i 0 , 2 i0-1 k 0 ). Let us consider a pyrel a 0 (i 0 , k 0 ). Its wavelet transform is:

w(i, k) = a 0 v(i, i 0 , k, k 0 )) (9) v(i, i 0 , k, k 0 ) is the pyramidal wavelet transform of h(i 0 , m + 2 i0-1 k 0 ).
Due to the decimation from a scale to the following one, this function depends on k -2 i0-i k 0 , up to scale i 0 . The pyrel amplitude is first estimated by:

a(i 0 , k 0 ) = w(i 0 , k 0 ) v(i 0 , i 0 , 0, 0) (10)
In the case of isolated pyrels, the amplitude is correctly estimated by the relation [START_REF] Tropp | Greed is good: algorithmic results for sparse approximation[END_REF]. If two pyrels at the same scale are too close, it is not possible to separate them directly by scanning the extrema of the wavelet transform. But, if the pyrels are not at the same scale, the largest scale pyrel can be considered as background for the smallest one, and a separation is then possible. The matching pursuit algorithm can proceed at each step on the extrema set. All the extrema having a SNR greater than the chosen threshold are simultaneously considered. The extrema which correspond to the extrema of the wavelet coefficients along the scale are kept, are called suprema. Thus, the pyrel identification consists in two steps:

1) Identification, scale by scale, of the local extrema of the wavelet transform which have a SNR greater than the threshold; 2) Identification of the suprema, corresponding to the extremum along the scale of the wavelet transform.

D. The pyrel coupling.

The pyrels are previously assumed to be separated. A pyrel is identified, its amplitude evaluated and the pyrel subtracted, and so one. But the pyrels are not orthogonal patterns. They are identified from the wavelet coefficients, and the wavelet functions are not also orthogonal. So a coupling exists between the pyrels. The signal is restored by the relation

F (k) = l a(i l , k l )h(i l , k -2 i l -1 k l ) ( 11
)
where l is the pyrel index, (l ∈ (1, L). The wavelet transform of the reconstructed signal is:

w(i, k) = l a(i l , k l )v(i, i l , k, k l ) (12)
The pyrel amplitudes being determined from their wavelet coefficients at the same location, both in scale and position, the following coupling equation results:

w(i l , k l ) = l ′ a(i l ′ , k l ′ )v(i l , i l ′ , k l , k l ′ ) (13)
An implicit procedure is done for inverting Equation [START_REF] Strang | Wavelets and Dilation Equations: A Brief Introduction[END_REF]. It takes into account the wavelet transform of the image reconstructed from the pyrels. A Van Cittert correction [START_REF] Van Cittert | Zum einfluss der spaltbreite auf die intensitäts verteilung in spektrallinen II[END_REF] of the amplitude is done by adding iteratively on each amplitude a term proportional to the difference between the original wavelet coefficient and its restored one.

E. The greedy algorithm from a pyramidal representation.

The pyrel representation from the pyramidal wavelet transform is derived from the previous analysis:

1) Computation of the pyramidal wavelet transform 2) Determination of the extrema of the wavelet transform at each scale 3) Determination of the suprema. 4) Estimation of the maximum SNR, Q max . A threshold equal to αQ max is chosen for the suprema selection. If Q max is little than a given threshold, the loop is stopped and the algorithm goes to step [START_REF] Chen | Atomic Decomposition by Basis Pursuit[END_REF]. 5) The detected wavelet coefficients w(i j , k j ), j ∈ (1, J), are multiplied by a factor which takes into account the ratio between the pyrel amplitude and the wavelet one. 6) The identified pyrels are added to the previous detected ones, which leads to a set {a (0) (i l , k l )} where l corresponds to the set of identified pyrels for all the iterations. 7) An inverse algorithm is done in order to determine the new set {a (∞) (i l , k l )} such that the signal reconstructed with these pyrels carries out after a wavelet transform which has the wavelet coefficients w(i l , k l ) at the pyrel locations. 8) An image is built from the pyrels. Its pyramidal wavelet transform is computed and it is subtracted to the original one. A test is done if significant wavelet coefficients still exist. If yes, the algorithm goes to step [START_REF] Usevitch | A tutorial on modern lossy wavelet image compression: foundations of JPEG[END_REF]. 9) If the restored signal has to be close to the original one, a baseline is added. The baseline is estimated as the difference between the last approximation of the signal and the last approximation of the pyrel representation.

As the sampling is progressively reduced, the baseline is reinterpoled by a set of correlations with 2h(n), after the insertion at each step a 0 between two values. At the largest scale, the baseline can be reduced to a constant.

F. The algorithm tuning.

The algorithm depends on different parameters. The low-pass filter. The choice of the low-pass filter is governed by the following constraints:

• Its length must be as short as possible in order to reduce the computations and the boundary artifacts.

• The corresponding scaling function φ(x) must be always positive and must display only one maximum, in order to associate a pyrel to a maximum.

• The filter must be symmetric, in order to get available pyrel positions.

• φ(x) must be as regular as possible, taking into account its length. This will avoid to introduce pyrels due to the irregularities of the scaling function. The binomial filter {h(n) = C n+l 2l }, for -l ≤ n ≤ l, corresponds to the centered B-spline of order 2l -1 for the scaling function. It fully satisfies the constraints. In our application, we use the filter for l = 2 that corresponds to the cubic centered B-spline. The α factor. In the algorithm description, an α factor has been introduced in order to select the suprema which have to be processed in parallel. In the program, it was set to 0.5. With this value it is assumed not to select suprema whose the values would be too contaminated by the new identified pyrels. If this factor decreases, the algorithm may consider the bumps of the coupling function. Experimentally, this factor appears also as a good compromise between sparsity (which increases with α) and computing time (which also increases with α). The break parameter in the inversion. The pyrel amplitudes are determined by inversion which takes into account the coupling matrix. This inversion is done by a Van Cittert iterative algorithm. The iterations are stopped when the highest residual is greater than rQ max . A value of r = 0.25 was experimentally a good compromise . The Van Cittert convergence factor. The Van Cittert iterative algorithm was designed in order to solve Y = AX, where A is smoothing operator. The iterations are written as:

X (n+1) = X (n+1) + β[Y -AX (n) ]. (14)
β is not necessary equal to 1. The stability needs 0 < β < 2, generally β = 1. The algorithm instability in case of a bad matrix conditioning is well known in deconvolution problems. After experiments, we set in the algorithm β = 0.5.

G. The monosign pyrel decomposition.

In the previous algorithm the pyrel amplitude could be either positive or negative, according to the supremum sign. Our scientific purpose was to build objects from pyrels. The more often astronomical images display positive objects drowned in a background. It is thus convenient, and easy, to modify the algorithm in order to force the detection of only positive (or only negative) peaks.

We note that the wavelet transform of a real negative peak leads to two positive peaks due to the bumps of the wavelet function. Therefore, in case of a positive decomposition, a peak which significantly corresponds to a negative object can be identified as positive structures. Nevertheless, the bump amplitudes are lower than the central peak and some significant negative peaks can remain after the decomposition. Tests on the residual wavelet coefficients can alarm that significant negative (or positive) peaks still remain.

H. Some properties.

The effect of translations. Due to the decimation, the decomposition is not shift-covariant. It is possible to get this covariance if no decimation is done from one scale to the following one. The algorithm complexity is increased. Dyadic dilation or contraction. The pyrels are generated by dilation of a generic pattern. If the signal is dilated by the same dilation, as a first approximation, the same pyrels would be identified at the above scales. But, in details some small changes may appear. For the contraction, the pyrel decomposition can differ sensibly from a simple shift along the scales. Non linearity of the decomposition. The algorithm is covariant with the scalar multiplication if the threshold is multiplied by the same scalar. Let us consider now two signals {F 1 (k)} and {F 2 (k)} leading to pyrel decompositions {a 1l (i l , k l )} and {a 2l (i l , k l )}. The signals are added:

• For pyrels in the signals at the same location, the resulting pyrel amplitude is the sum of their amplitudes. If they have not the same sign, the new amplitude may be less than the threshold and the pyrel will not be detected.

• Added pyrels may merge into a larger one, more significant.

• If pyrels of each signal have compact separate support, the adding is kept. If not, the decomposition is modified by the pyrel coupling. The pyrel decomposition of the result is thus generally different of the union of the pyrels coming from each decomposition. Due to the thresholding and the coupling between the pyrels, this decomposition is not linear. Invariance to a baseline addition. This algorithm is built to be invariant to the addition of a constant. More generally, the algorithm may be invariant to the adding of a non constant baseline whereas its wavelet transform does not modify the detection at the chosen threshold, in the scale range. One pyrel identification. Let us consider a signal composed only of a pyrel a 0 (i 0 , k 0 ). If the amplitude is sufficiently high, the algorithm only detects it. The other extrema (due to the wavelet bumps) are removed taking into account the suprema identification rules. The wavelet amplitude is converted into one pyrel. The pyrel is computed and subtracted. The residuals are null, apart the computational errors. A single pyrel is correctly restored. Multiple pyrels recognition. Let us consider now a signal composed by L pyrels {a l (i l , k l )}. The question is to know if the algorithm would restore these pyrels. Evidently this could be possible only if the wavelet coefficients related to the pyrels at the same scales and locations would be greater than the threshold.

But even if this condition is satisfied, there are reasons for which the full identification is not possible:

• The merging of pyrels into larger ones;

• The coupling between them which leads to increase the number of identified pyrels. Nevertheless, it is clear that the recognition of input pyrels works well their supports are disjoined. The decomposition algorithm can also recognize pyrels even if this condition is not satisfied. The main difficulty lies into its capability to detect the pyrels at their correct location (position and scale) during the identification phase.

A basic condition to recognize a set of pyrels is their linear independence. If the pyrels are not independent, their Gram-Schmidt matrix is singular and the inversion fails. Stability of the decomposition. Let us consider a signal which leads to a set of pyrels {a

(0) l (i (0) l , k (0)
l )} after decomposition with a threshold T . The signal is reconstructed from the pyrels and the decomposition algorithm is applied on it, with the same threshold. The decomposition algorithm is stable on this signal if the restored set {a

l (i (1) l , k (1)
l )} is equal to the initial one.

An exact stability is not generally strictly reached, taking into account the computational accuracy. The decomposition is also sensitive to the threshold. In the case of pyrels having amplitudes close to the threshold, computational errors may also lead to remove these features. Numerical experiments showed that stability was not assumed on studied signals. The restored signals are quite identical, but the number of pyrels, their positions, their scales and their amplitudes may sensibly differ.

I. Relation to a variational approach.

This instability shows that this greedy algorithm does not bring the global minimum of a given functional. The algorithm carries out a decomposition interesting for a given analysis. An associated variational approach, if it exists, would improve its stability. The data attachment. This attachment is generally done by the l 2 distance between the observed data and the restored ones. Here, the removal of the background avoids the application of this principle. The algorithm identifies pyrels from the extrema in the wavelet space, up to a given scale. The identification is based on a thresholding. The data attachment is thus connected to the pseudo-norm:

W ∞,I = max i=1,I,k=1,K | w(i, k) -w(i, k) σ i | (15)
The condition W ∞,I ≤ t indicates that a variational version would be connected to a minimax approximation [START_REF] Powell | Approximation theory and methods[END_REF] in the wavelet transform space (WTS). A function in this space w(i, k) is decomposed with the set {v(i, i p , k, k p )} , where p covers all the scales and positions for a pyrel. The prior condition. Today the prior on the decomposition is often linked to the ℓ 1 norm, i.e. the sum of the absolute coefficient values. The basis pursuit corresponds to this variational constraint [START_REF] Chen | Atomic Decomposition by Basis Pursuit[END_REF]. An algorithm which minimizes ℓ 1 taking into account W ∞,I would converge to a sparse decomposition. Nevertheless, the proposed algorithm is based on a direct identification of the pyrels, starting from the most significant ones. The goal is to reduce the number of elements allowing one to represent the image such that W ∞,I reaches a given value. So, the prior is not connected to the sum of the absolute values of the coefficients, but to their number.

It was shown previously that the results were not stable, but the variations correspond to quite insignificant coefficients resulting of progressive approximations. An algorithm which constrains to minimize the number of pyrels with W ∞,I ≤ t, t being a given threshold, would bring the stability. Greedy algorithms which proceed by decreasing the identification thresholds are a natural way to solve the problem, but they do not allow to carry out the stable solution which would correspond to the global minimum in case of non convexity.

Minimax approximations are extensively used in order to approximate functions by rational approximations. The exchange algorithm allows one to get the decomposition, if the Haar condition on the decomposition is satisfied [START_REF] Powell | Approximation theory and methods[END_REF]. In the present variational problem, this condition is not satisfied, so that another class of algorithm has to be developed.

J. Two dimensional algorithm.

The 2D multiresolution pyramid. The multiresolution pyramid is defined by:

F (i, k) = 1 4 i < f (x), φ( x 2 i -k, y 2 i -l) > . (16)
The following recursive relation is derived in case of separate variables:

F (i + 1, k, l) = n h(n)h(m)F (i, 2k + n, 2l + m). ( 17
)
The derived patterns. Taking into account the separation between the variables, the coefficients can be also written as:

F (i + 1, k) = p=p1,p2 q=q1,q2 h(i, p)h(i, q) F (2 i+1 k + p, 2 i+1 k + q), (18)
the filters h(i, p) being the 1D pyrels.

The identification from the wavelet transform. Similarly, the wavelet transform is built on the difference between two successive approximations. The suprema of the resulting transform are identified, taking into account a threshold, derived from the noise level.

In order to avoid identifying first spurious pyrels, the threshold is progressively decreased, with the same rule that the one used for the 1D. The algorithm. Taking into account the previous considerations, the 2D algorithm is simply copied from the 1D one. The image can be decomposed with only positive (or negative) pyrels, taking into account the suprema sign.

III. AN APPLICATION ON AN ASTRONOMICAL IMAGE.

On Figure 1 left, the image of the planetary nebula NGC40, taken in the near infrared, is plotted. This image is characteristic of the astrophysical images. A bright central star illuminates a shell, which splitted in two parts. The image contains star-like objects and extended diffuse ones. On Figure 2 the image of the planetary nebula NGC40 after the application of this decomposition algorithm, with the unsigned analysis. Here, the number of pyrels corresponds to 3.1% of the number of pixels.

On Figure 3 the image is reconstructed from the decomposition based on only positive pyrels. Now, the number of pyrels corresponds to only 2.6% of the number of pixels.

Using this decomposition, the reduction of the number of significant coefficients is a factor 2 compared to the thresholded orthogonal discrete wavelet transform. Not only the sparsity is seriously increased, but also the reconstruction does not display the wavy artifacts connected to the inverse wavelet transform.

IV. CONCLUSION

In the present paper, we propose a new way to decompose a signal or an image. Its main specificity is the local background

Pyramidal postive decomposition Residuals

Fig. 3. The reconstructed image of the planetary nebula NGC40 from positive pyrels identified with the pyramidal algorithm (left). The residual image is quite a white noisy image (right). Nevertheless in the central part quite significant structural features can be discerned.

removal. In the framework of a multiresolution analysis that leads to examine the suprema of the wavelet transform in order to reconstruct with scalet patterns, here called pyrels. The matching pursuit is then done simultaneously for a set of coefficients, while their amplitudes are greater than a threshold which progressively decreases up to a given level. Sufficiently separated suprema are kept at each step. An amplitude correction is also made in order to recover the observed wavelet coefficients at the pyrel locations.

The algorithm (MSMPPy1) was presented for a 1D signal, with the use of an pyramidal wavelet transform. The application of the undecimated wavelet transform increased the complexity, but it leads to a shift invariant transform (MSMPAT1). The two-dimensional algorithms (MSMPPy2 and MSMPAT2) are their natural extensions to the two-dimensional field.

MSMPAT1 is quite fast and it can be applied on a large size signal. Nevertheless, MSMPPy1 is more appropriate to process long series. For an image, MSMPPy2 is the useful tool for the analysis of large images, but MSMPAT2 brings a better description. A careful analysis of small images is more convenient with this tool. For these four algorithms a signed decomposition (positive or negative) can be done.

The pyramidal algorithms carry out a sparse decomposition. For example, on the planetary nebula image, a gain around 2 was obtained in the number of coefficients, compared to a classical wavelet thresholding. The different experiments on astronomical images shows that this gain was the more often higher. The image compression is thus a direct application of these algorithms.

These algorithms were built for the analysis of multiband astrophysical images. For each band, it is easy to determine the amplitude of each pyrel taking into account the same identification set. The merging of pyrels into objects allows one to give a full description of the images. A complete multiband vision was then derived [START_REF] Bijaoui | Fusion images and Analysis of multiband astronomical images[END_REF].

Fig. 1 .

 1 Fig. 1. The infrared image of the planetary nebula NGC40.

Fig. 2 .

 2 Fig.2. The reconstructed image of the planetary nebula NGC40 from pyrels identified with the pyramidal algorithm (left). The residual image looks like to a white noisy image (right).


Cristiano Maria Verrelli

Riccardo Marino

Patrizio Tomei

Gilney Damm

Nonlinear Robust Coordinated PSS-AVR Control

Keywords: Nonlinear control, Synchronous generators, Transient stability, Voltage regulation, Robust coordinated PSS-AVR

come   L'archive ouverte pluridisciplinaire

I. INTRODUCTION

The problem of power systems stability [START_REF] Anderson | Power System Control and Stability[END_REF], [START_REF] Sauer | Power System Dynamics and Stability[END_REF] is still widely studied, according to the latest related theoretical results in [START_REF] Barabanov | Conditions for almost global attractivity of a synchronous generator connected to an infinite bus[END_REF] and [START_REF] Weitenberg | Robustness to noise of distributed averaging integral controllers in power networks[END_REF], [START_REF] Weitenberg | Exponential convergence under distributed averaging integral frequency control[END_REF], [START_REF] Weitenberg | Robust decentralized frequency control: a leaky integrator approach[END_REF], [START_REF] Weitenberg | Robust decentralized secondary frequency control in power systems: merits and trade offs[END_REF]. Transient stability is to be enhanced in the presence of perturbations such as changes in load, outages of powers plants, failures in transformer substation and power lines. Electromechanical oscillationsnamely, local, inter-machine, and inter-area oscillations -have to be damped, while terminal voltage is to be regulated to its reference value. Linear and local controllers are currently employed. In practice, Linear Automatic Voltage Regulators (AVRs) are used in conjunction with Power System Stabilizers (PSSs) (see [START_REF] Bourlés | Analysis and design of a robust coordinated AVR/PSS[END_REF], [START_REF] Dudgeon | The effective role of AVR and PSS in power systems: frequency response analysis[END_REF], [START_REF] Dyśko | Enhanced power system stability by coordinated PSS design[END_REF], [START_REF] Law | Robust co-ordinated AVR-PSS design[END_REF], [START_REF] Machowski | Decentralized stability-enhancing control of synchronous generator[END_REF], [START_REF] Quinot | Robust coordinated AVR+PSS for damping large scale power systems[END_REF]). Those controllers are typically numerically designed and tuned on the basis of the classical reduced model of a synchronous generator connected to an infinite bus, whose dynamics are linearized around the (unknown) desired constant operating condition that guarantees terminal voltage regulation. Such a linear model will be referred, throughout this technical note, to C.M. Verrelli, R. Marino, P. Tomei are with the Electronic Engineering Department of the University of Rome Tor Vergata, Via del Politecnico 1, 00133 Rome, Italy, E-mail (of the corresponding author): verrelli@ing.uniroma2.it. G. Damm is with COSYS-LISIS, Univ Gustave Eiffel, IFSTTAR, F-77447 Marne-la-Vallée, France.

as LM. Even though satisfactory results are achieved by the "Desensitized Four Loops Regulator" DFLR in [START_REF] Bourlés | Analysis and design of a robust coordinated AVR/PSS[END_REF] and [START_REF] Quinot | Robust coordinated AVR+PSS for damping large scale power systems[END_REF] (here referred to as robust coordinated PSS-AVR in the whole paper), two open problems remain: P 1 : in order to obtain satisfactory performance over a wide range of operating conditions, different behaviours of nonlinear power systems in different operating conditions require different control objectives and consequently different linear control gain scheduling; P 2 : the absence of nonlinear feedback terms in non-switching controllers may prevent the enlargement of the machine stability region. Two different strategies are then in order. Switching Linear Controllers: the first strategy consists in using different linearization-based controllers with switching actions depending on the different operating conditions (see [START_REF] Guo | Global transient stability and voltage regulation for power systems[END_REF] and references therein for applications of gain scheduling, fuzzy control and multi-model techniques). Nonlinear Non-Switching Controllers: the second strategy consists in resorting to nonlinear non-switching controllers that incorporate nonlinear feedback terms, simultaneously allowing for the automatic gain scheduling and enlargement of the machine stability region. The structures of such controllers however are rather complex and require additional measurements and parameters. Examples are provided: i) a state feedback solution is proposed in [START_REF] Damm | Transient stabilization and voltage regulation of a synchronous generator[END_REF] (see also [START_REF] Marino | Robust adaptive transient stabilization of a synchronous generator with parameter uncertainty[END_REF], [START_REF] Marino | A global state feedback output regulating control for uncertain systems in strict feedback form[END_REF]), in the presence of uncertainties in almost all system parameters; ii) output feedback solutions are reported in [START_REF] Verrelli | Output feedback transient stabilization and voltage regulation of a synchronous generator[END_REF] and [START_REF] Damm | Adaptive nonlinear output feedback for transient stabilization and voltage regulation of power generators with unknown parameters[END_REF], in the presence of transmission line impedance and mechanical power uncertainties. In contrast to the robust coordinated PSS-AVR strategy, the strategy adopted in the aforementioned papers (see also [START_REF] Jiao | Adaptive controller design for a synchronous generator with unknown perturbation in mechanical power[END_REF]) does not consider the terminal voltage as the natural controlled output. Rather, it artificially transforms the power angle into the main controlled output. As a consequence, the (available) terminal voltage regulation error is not directly used in the feedback action, with either the power angle being consequently assumed to be directly measured or the knowledge of critical parameters being required. This technical note presents two contributions. The first contribution provides new interpretations for the "Desensitized Four Loops Regulator" of [START_REF] Bourlés | Analysis and design of a robust coordinated AVR/PSS[END_REF], by viewing it as a robust stabilizer for the stably invertible system LM with respect to a suitably chosen output function. The second contribution takes advantage from the previously defined output and presents a new nonlinear controller that generalizes the robust coordinated PSS-AVR -namely, its linear action coincides with the one provided by the robust coordinated PSS-AVR -, with, additionally, no use of the mechanical input power. In fact, such a nonlinear generalization just endows the robust coordinated PSS-AVR with: i) a nonlinear term that relies on the power angle (being computed as in [START_REF] Verrelli | Output feedback transient stabilization and voltage regulation of a synchronous generator[END_REF] or being measured via now available wide-area measurement technologies), allowing for an automatic gain scheduling; ii) nonlinear robustifying terms on the back-stepping-based electrical power tracking error, with the aim -whose effectiveness is illustrated by a numerical analysis -of enlarging the machine stability region. The simplicity of control structure and tuning procedure of the linear design is maintained: just one integrator is involved, while all model parameters are uncertain and are not estimated by the controller. From a theoretical point of view, such a nonlinear contribution simultaneously improves: i) the results in [START_REF] Gao | A nonlinear control design for power systems[END_REF], [START_REF] Guo | Global transient stability and voltage regulation for power systems[END_REF], [START_REF] Wang | Robust decentralized nonlinear controller design for multimachine power systems[END_REF], [START_REF] Zhu | A new nonlinear voltage controller for power systems[END_REF], since transient stabilization and voltage regulation tasks are simultaneously accomplished, with no use of robust design strategies that require the existence of suitable stabilizing matrices and with no need of model-knowledgebased feedback linearizing precompensation-actions or switching between transient stabilizers and voltage regulators; ii) the output feedback results in [START_REF] Verrelli | Output feedback transient stabilization and voltage regulation of a synchronous generator[END_REF] (see also [START_REF] Damm | Adaptive nonlinear output feedback for transient stabilization and voltage regulation of power generators with unknown parameters[END_REF] and references therein), since, even when the power angle profile is computed as in [START_REF] Verrelli | Output feedback transient stabilization and voltage regulation of a synchronous generator[END_REF], only the values for the generator direct axis reactance and the voltage of the infinite bus are here required to be available, whereas several model parameters (including the damping coefficient) are assumed to be known in [START_REF] Verrelli | Output feedback transient stabilization and voltage regulation of a synchronous generator[END_REF]; iii) the state feedback results in [START_REF] Damm | Transient stabilization and voltage regulation of a synchronous generator[END_REF] (and all references therein), since several parameter estimates are not required, thus reducing the control complexity. Realistic simulations concerning a benchmark scenario show that the nonlinear controller preserves the stability when a short circuit occurs, whereas instability arises when the standard PSS-AVR is used.

II. PROBLEM FORMULATION AND DFLR CONTROL

In this section we present the generator dynamic model, formulate the problem and introduce the robust coordinated PSS-AVR control (DFLR control).

A. Dynamic model

The third order dynamic model of a synchronous generator connected to an infinite bus (see [START_REF] Damm | Adaptive nonlinear output feedback for transient stabilization and voltage regulation of power generators with unknown parameters[END_REF] for its description and [START_REF] Verrelli | Robust transient stabilization problem for a synchronous generator in a power network[END_REF], [START_REF] Barabanov | Conditions for almost global attractivity of a synchronous generator connected to an infinite bus[END_REF] for more general models) is given by:

δ = ω ω = - D H ω - ω s H P e + ω s H P m (1)
Ṗe = -θ 1 P e + θ 2 ω sin 2 (δ) + P e ωcotg(δ)

+ θ 3 sin(δ)u f V t = θ 4 P 2 e sin 2 (δ) + θ 5 + θ 6 P e cotg(δ) 1 2
, in which: δ (rad) denotes the generator power angle with respect to the infinite bus rotating at synchronous speed ω s ; ω = ω g -ω s (rad/s) is the difference between the generator angular speed ω g and the synchronous speed ω s ; P e (p.u.) denotes the active electrical power delivered by the generator to the infinite bus; u f (p.u.) is the input to the Silicon Controlled Rectifier amplifier of the generator; P m (p.u.) is the mechanical input power, V t (p.u.) is the generator terminal voltage. The reader is referred to [START_REF] Damm | Adaptive nonlinear output feedback for transient stabilization and voltage regulation of power generators with unknown parameters[END_REF] for the physical meaning of D, H and the six positive parameters θ i > 0, 1 ≤ i ≤ 6.

Model (1) turns out to be the well-known single-machineinfinite-bus representation, which, though of low dynamic order, well captures the main (small-signal) characteristics of a synchronous generator in a power network in the presence of parameter variations (see [START_REF] Bourlés | Analysis and design of a robust coordinated AVR/PSS[END_REF], [START_REF] Quinot | Robust coordinated AVR+PSS for damping large scale power systems[END_REF]). Now, in order to avoid the physical singularities at δ = 0 and δ = π, the operation of the system is restricted to an arbitrarily large compact (connected) set B contained in the open set: D = (δ, ω, P e ) : (δ, ω, P e ) ∈ (0, π) × R × R + .

B. Problem formulation

Definition 1 (transient stabilization & voltage regulation): Consider the synchronous generator described by the dynamic model [START_REF] Anderson | Power System Control and Stability[END_REF], with its operation being restricted to B. Assume that all the model parameters are constant and uncertain. The (robust) transient stabilization and voltage regulation problem for the synchronous generator consists in exponentially regulating to zero: i) the relative speed ω (transient stabilization); ii) the error between terminal voltage V t and its reference value V tr = 1 (p.u.) (voltage regulation).

The problem formulation of Definition 1 is well-posed, since the desired operating condition (at which voltage regulation is achieved) exists. It is given by δ = δ s , ω = 0, P e = P m , where (see [START_REF] Damm | Adaptive nonlinear output feedback for transient stabilization and voltage regulation of power generators with unknown parameters[END_REF] for the physical meaning of the parameters):

δ s = arccotg - V s 2 X d X s P m (X s + X d ) + V s 2 V tr 2 X s 2 P m 2 -1
belongs to the open set (0, π/2) and guarantees terminal voltage regulation, namely:

θ 4 P 2 m sin 2 (δ s ) + θ 5 + θ 6 P m cotg(δ s ) 1 2 = V tr . (2)
It corresponds to a stable equilibrium for system [START_REF] Anderson | Power System Control and Stability[END_REF] fed by the open loop control u f s = θ 1 P m θ -1 3 sin(δ s ) -1 . According to [START_REF] Damm | Transient stabilization and voltage regulation of a synchronous generator[END_REF], there is also another value δ u ∈ π/2, π for δ (see also [START_REF] Barabanov | Conditions for almost global attractivity of a synchronous generator connected to an infinite bus[END_REF]) that guarantees terminal voltage regulation. It however corresponds to an unstable equilibrium for system [START_REF] Anderson | Power System Control and Stability[END_REF] fed by the open loop control

u f u = θ 1 P m θ -1 3 sin(δ u ) -1 .

C. Robust coordinated PSS-AVR control (DFLR)

Let us introduce the following regulation errors: ω = ω, Pe = P e -P m , Ṽt = V t -V tr . The classical Desensitized Four Loops Regulator (DFLR) (or, equivalently, the robust coordinated PSS-AVR) in [START_REF] Bourlés | Analysis and design of a robust coordinated AVR/PSS[END_REF] [see also (2)-(3) in [START_REF] Quinot | Robust coordinated AVR+PSS for damping large scale power systems[END_REF]] can be then written as:

u f (t) = k ω ω(t) -k p Pe (t) -k v Ṽt (t) -k I t 0 Ṽt (τ )dτ, (3)
in terms of suitable control parameters k ω , k p , k v , k I . The above DFLR is apparently constituted by the combination of: i) the (transient) stabilizing control action

u f,ts (t) = k ω ω(t) -k p Pe (t) (4)
that is represented, according to the second equation of ( 1), by the transfer function (acting on the Laplace transform of the relative speed ω):

W ts (s) = k ω + k p Hω s -1 s + DH -1 ; (5)
ii) the voltage regulating control action

u f,vr (t) = -k v Ṽt (t) -k I t 0 Ṽt (τ )dτ (6)
that is represented by the transfer function (acting on the Laplace transform of the terminal voltage regulation error Ṽt ):

W vr (s) = -(k v + k I /s) . (7)
As shown by [START_REF] Bourlés | Analysis and design of a robust coordinated AVR/PSS[END_REF], such DFLR can be put, by using straightforward approximations that do not significantly deteriorate the performance of the controller, in the standard AVR-PSS structure IEEE ST1A + PSS1A reported in Figures 2 and4 of [START_REF] Bourlés | Analysis and design of a robust coordinated AVR/PSS[END_REF], with ( 5) and ( 7) mimicking the actions of the PSS and AVR, respectively.

III. NEW INTERPRETATIONS FOR THE DFLR

Define, besides the already defined regulation errors (ω = ω, Pe = P e -P m , Ṽt = V t -V tr ), the power angle regulation error: δ = δ-δ s . Write the first order approximated expression of Ṽt (around the desired operating condition [δ s , 0, P m ] T ) as Ṽt = -β 0 + β 1 Pe , where β 0 and β 1 turn out to be positive over the generator operating region (see [START_REF] Zhu | A new nonlinear voltage controller for power systems[END_REF]). We preliminarily state Lemma 1, whose results will be used in the reminder of this technical note.

Lemma 1: The linearized dynamics LM of the synchronous generator around the desired operating condition [δ s , 0, P m ] T is of relative degree one and minimum phase with respect to the output

y = Ṽt -µω + ν Pe , (8)
where µ and ν are reals satisfying

ω s µ + (β 1 + ν)D > 0 & β 1 + ν > 0.
Proof: Starting from (1), consider the linearized dynamics LM in terms of the state variables: δ, ω, Pe and the output y. By mimicking the first steps in [START_REF] Verrelli | Output feedback transient stabilization and voltage regulation of a synchronous generator[END_REF], perform the change of coordinates:

z 1 = δ, z 2 = ω, z 3 = -DH -1 ω -ω s H -1
Pe and thus obtain the system expressed in the new coordinates. It is of relative degree one and in reachability form: its transfer function (between the output y and the input u ω s θ 3 H -1 sin(δ s )(u f -u f s )) can be derived to be:

W M (s) = (β1+ν)H ωs s 2 + µ + (β1+ν)D ωs s + β 0 s 3 + D H -α 2 s 2 + ωs H α 1 -D H α 2 s + ωs H α 0 where α 0 , α 1 , α 2 are the reals: α 0 = θ 3 u f s cos(δ s ), α 1 = θ 2 sin 2 (δ s ) + P m cotg(δ s ), α 2 = -θ 1 depending
on the generator operating condition. Owing to the positive nature of the coefficients characterizing the numerator of W M (s), the minimum phase property is guaranteed.

We now state Theorem 1.

Theorem 1: Consider the linearized dynamics LM of the synchronous generator around the operating condition [δ s , 0, P m ] T . Let µ, ν satisfy the constraints of Lemma 1. Then the application of the (µ, ν)-DFLR (k y , k yI > 0):

u f,µ,ν (t) = -k y y(t) -k yI t 0 y(τ )dτ (9)
achieves transient stabilization and voltage regulation, provided that k y is sufficiently large. Proof: Equivalently express the generalized (µ, ν)-DFLR (9) as:

u f,µ,ν (t) = -k y y(t) + ûfs (t) ufs (t) = -k yI y(t), ûfs (0) = 0. ( 10
)
Consider the transfer function W M (s). For the sake of compactness: i) denote by b 1 , b 2 , b 3 the positive coefficients of its numerator and let them constitute the components of a column vector b, with the zeroes of the polynomial p(s) = b 1 s 2 + b 2 s + b 3 all belonging to C -; ii) denote by a 1 , a 2 , a 3 the coefficients of its monic denominator and let them constitute the components of the column vector a. Accordingly take the realization of W M (s) in observability form:

ẋ =   0 1 0 0 0 1 0 0 0   x + ay + bu A c x + ay + bu y = [1, 0, 0]x. (11)
Recalling that b 1 > 0, the linear change of coordinates

[η T , y] [η 1 , η 2 , y] -with η i = x i+1 -b i+1 b 1 -1 x 1 , 1 ≤ i ≤ 2, y = x 1 -
can be performed and stability arguments similar to those used in [START_REF] Marino | Learning control for nonlinear systems in output feedback form[END_REF], once specialized to the case of constant disturbance compensation, can be used to prove the thesis.

The following Corollary demonstrates that the above (µ, ν)-DFLR ( 9) is related to a reparameterization of the DFLR (3).

Corollary 1: Consider the linearized dynamics LM of the synchronous generator around the operating condition [δ s , 0, P m ] T . Let Pe0 , ω0 , Ṽt0 denote the initial conditions for Pe (t), ω(t), Ṽt (t) at time t = 0. Define the constant:

σ 0 = k yI µ + νD ω s β 1 β 0 Pe0 - Ṽt0 β 0 + k yI νH ω s ω0 . ( 12
)
Then, the (µ, ν)-DFLR [START_REF] Guo | Global transient stability and voltage regulation for power systems[END_REF] with ûfs (0) = σ 0 becomes the DFLR (3), with k v , k ω , k p satisfying:

k ω = k y µ + k yI ω -1 s Hν, k p = k y ν -k yI β 1 β -1 0 µ + νDω s -1 , k v = k y + k yI β -1 0 µ + νDω s -1 , k I = k yI . Proof:
Use the δ and ω dynamics for the computation of the y-integral and write δ as β -1 0 (β 1 Pe -Ṽt ). Once Lemma 1 and Corollary 1 have been established, a direct guideline in the control parameters tuning for (3) can be easily obtained: it involves just one parameter to be set at last, in accordance with the following corollary.

Corollary 2 (limit behaviour): Consider the linearized dynamics LM of the synchronous generator around the operating condition [δ s , 0, P m ] T . Then, the application of the (µ, ν)-DFLR (9) with k yI /k y being constantly set to z, for k y tending to +∞, makes: i) two of the closed loop eigenvalues converge to the two open loop zeros of W M (s) (not explicitly depending on X d , T d0 , K c ); ii) one of the closed loop eigenvalues converge to -z; iii) one of the closed loop eigenvalues tend to -∞ (along the real axis).

Proof: Consider LM, expressed, in accordance with the proof of Lemma 1, in terms of the coordinates z 1 , z 2 , z 3 and rewrite u f,µ,ν in terms of such coordinates. Define the new variable

Hz 4 = -ω s α 0 z 1 -(ω s α 1 -Dα 2 ) z 2 -(D -Hα 2 ) z 3 -ω s θ 3 sin(δ s )(u f -u f s ) (13)
and compute its dynamics. The proof relies on the root locus analysis and on the structure of the corresponding characteristic polynomial for the error system in z 1 , z 2 , z 3 , z 4coordinates, which reads

q z (s) = m 1 (s) + k y ω s sin(δ s )θ 3 H -1 m 2 (s), (14)
where the roots of the polynomials m 1 (s) and m 2 (s) denote the open loop poles and zeroes, respectively (recall W M (s) and the Laplace transform of u f,µ,ν in ( 9)). Such open loop zeroes, in accordance with the definition of β 0 , β 1 and θ 4 , θ 6 , δ s , do not explicitly depend on X d , T d0 , K c .

The proof of such a corollary clarifies that non-zero µ and ν in ( 9) are assigned to guarantee satisfactory transient performance, and in particular, to satisfactorily damp down speed oscillations: if µ = ν = 0, then, for arbitrarily large k y (and k yI /k y being constantly set to z), two closed loop eigenvalues are attracted by the open loop zeroes of W M (s), which, for such µ = ν = 0, guarantee a non-satisfactory damping since D is typically very small. On the other hand, Corollary 2 suggests to choose µ, ν, k yI /k y to robustly place the open loop zeroes in the left half complex plane C - (independently of X d , T d0 , K c ), and then to largely increase k y to achieve transient stabilization and voltage regulation 1 . The different parameterization of the DFLR (3) rather allows, to arbitrarily impose, once (3) is applied to the error system in z 1 , z 2 , z 3 , z 4 -coordinates, the four closed loop eigenvalues through the four control parameters k ω , k p , k v , k I (to be simultaneously chosen depending on the specific operating condition and all the system parameters).

IV. NONLINEAR GENERALIZATION OF THE DFLR

In this section, we consider the nonlinear model (1) in place of the linearized one LM. We accordingly present the design of a new nonlinear controller that generalizes, to the nonlinear model ( 1), the design of the robust coordinated PSS-AVR of Sections II-III, while additionally avoiding the use of the mechanical input power. The proposed nonlinear 1 It is interesting to notice that the industrial coordinated PSS-AVR of Figure 5 of [START_REF] Bourlés | Analysis and design of a robust coordinated AVR/PSS[END_REF] can be, in turn, obtained through the application of a slight variant of the approach underlying the proof of Theorem 1. Such a variant is based on Lemma 4 in [START_REF] Mareels | A simple selftuning controller for stably invertible systems[END_REF]. In particular, while the controller in Theorem 1 (described by a proper rational function) first feeds back the output error y to stabilize the (minimum phase) error system and then includes an integral action on y being responsible of generating a converging estimate of the uncertain constant u f s , such a different controller relying on Lemma 4 in [START_REF] Mareels | A simple selftuning controller for stably invertible systems[END_REF] (described by a strictly proper transfer function) first includes an integrator in its transfer function and then stabilizes the resulting error feedback through a proper rational function. In this respect, notice that Lemma 4 in [START_REF] Mareels | A simple selftuning controller for stably invertible systems[END_REF] is correct for d ≤ 4, whereas, for d ≥ 5, a typo concerning g 0 in place of q d appears (which the counterexample of [START_REF] Hoagg | Direct adaptive dynamic compensation for minimum phase systems with unknown relative degree[END_REF] to Lemma 4 of [START_REF] Mareels | A simple selftuning controller for stably invertible systems[END_REF] just relies on).

generalization of the robust coordinated PSS-AVR is given by: u f = -k p (P e -P er ) -k r Φ(δ, P e )(P e -P er ) + λ ρ sin(δ) Φ(δ, P e ) = sin2 (δ) 2 + [P e cotg(δ)]

2

P er = k pω ω + ρ (15) ρ = -µ v Ṽt + µ ω ω,
where: k p , k r , λ, k pω , µ v are suitable positive control parameters; µ ω is a suitable non-negative control parameter satisfying max{D, µ ω } > 0.

Before stating the main result of this section (namely, Theorem 2), the following comments are in order.

Comment A. The nonlinear robust coordinated PSS-AVR (15) involves the injection of the power angle to provide an automatic gain scheduling and to handle the severe disturbances and contingencies characterizing the operation of power systems. We can provide two interpretations for the above controller (15): i) the state feedback one, when the power angle is measured via now available wide-area measurement technologies (see [START_REF] Gonzalez-Torres | A novel distributed supplementary control of Multi-Terminal VSC-HVDC grids for rotor angle stability enhancement of AC/DC systems[END_REF]); ii) the output feedback one, when the values of δ and X s are computed through the (δ, X s )formulas derived in [START_REF] Damm | Adaptive nonlinear output feedback for transient stabilization and voltage regulation of power generators with unknown parameters[END_REF] (Q is the reactive power):

δ = arccotg - V s 2 X d X s P e (X s + X d ) ± V s 2 V t 2 X s 2 P e 2 -1 X s = -QV 2 s ± Q 2 V 4 s -(Q 2 + P 2 e )(V 2 s -V 2 t )V 2 s Q 2 + P 2 e
that just require the knowledge of the generator direct axis reactance and the voltage of the infinite bus 2 . With this respect, recall that X s may undergo sudden variations when electrical perturbations, such as faults on the transmission line, occur. Comment B. Set µ ω = 0 and ρ(0) = P m . Then, when higher order robustifying actions are neglected [w.r.t. the regulation and estimation errors δ, ω, Pe , P m -ρ], (15) reduces to ( kp = k p + k r Φ(δ s , P m )):

u f = 1 sin(δ) -kp Pe + kp k pω ω -λµ v Ṽt -kp µ v t 0 Ṽt (τ )dτ
that coincides with the DFLR (3):

u f = 1 sin(δ s ) -kp Pe + kp k pω ω -λµ v Ṽt -kp µ v t 0 Ṽt (τ )dτ
once δ appearing at the denominator is replaced by δ s and incorporated into the control gains. In other words, the nonlinear robust coordinated PSS-AVR (15) exhibits different linearizations around different operating conditions. Those different linear actions can be seen as equivalent actions of local robust coordinated PSS-AVR linear controllers with different control parameters. The controller ( 15) can be thus viewed as a nonlinear version of a simple combination between a robust coordinated PSS-AVR and an automatic smooth nonlinear gain scheduling procedure. The operating status of the controller is automatically changed when operating conditions vary, so that a smooth transfer between the actions of local robust coordinated PSS-AVR linear controllers is performed.

We are now able to present Theorem 2. Theorem 2: Consider the nonlinear robust coordinated PSS-AVR [START_REF] Mareels | A simple selftuning controller for stably invertible systems[END_REF]. Besides δ = δ -δ s and ω = ω, define the new error coordinates: e p = P e -P er -P ess , Pm = P m -P ess -ρ in terms of the suitable constant

P ess = -θ 1 P m θ 3 [k p + k r Φ(δ s , P m )] . ( 16
)
Then transient stabilization and voltage regulation is achieved.

In particular, the origin of the error system in the ( δ, ω, -(D+

ω s k pω )H -1 ω+ω s H -1
Pm , e p )-coordinates -corresponding to the desired operating condition [δ s , 0, P m ] T -, is exponentially stable, for sufficiently large k p and k r and for suitable µ ω , µ v , k pω robustly placing in C -the roots of the polynomial (not explicitly depending on

X d , T d0 , K c ): Π A (s) = s 3 + ((D + ω s k pω )H -1 + β 1 µ v )s 2 + (β 1 DH -1 µ v + ω s H -1 µ ω )s + ω s β 0 H -1 µ v .
Proof: First compute the δand the ωdynamics by defining λ ω = (D + ω s k pω )H -1 . According to the first part of Section III, write

Ṽt = -β 0 δ + β 1 Pe + h V ( δ, Pe ), (17)
where h V (•, •) collects the higher order terms in the regulation errors δ and Pe . By writing Pe = e p + k pω ω -Pm and by recalling the definition of P er , the Ṽt -expression, in the new error coordinates, becomes:

Ṽt = -β 0 δ + β 1 k pω ω -β 1 Pm + β 1 e p + η V ( δ, ω, Pm , e p )
where η V ( δ, ω, Pm , e p ) h V ( δ, Pe ) collects, by definition, the higher order terms in the new errors δ, ω, Pm , e p . The injection of the output y into the Pm -dynamics (no knowledge of the mechanical input power is here assumed, so that ν = 0 has to be taken in [START_REF] Gao | A nonlinear control design for power systems[END_REF]) is thus able to modify the last row of the matrix characterizing the ( δ, ω, Pm )error subsystem. Accordingly write the error dynamics in terms of the new error variables:

z 1 = δ, z 2 = ω, z 3 = -λ ω ω + ω s H -1
Pm . The aim is to recover the previously presented arguments concerning the error system in z-coordinates, once: i) the value of D has been artificially modified through the k pω -based prefeedback action (with the aim of obviating the absence of ν); ii) ρ has been endowed with the role of u f and y has been fed back with the suitable negative sign in order to stabilize, with the desired eigenvalues, the resulting third-order error subsystem. In this sense, the polynomial Π A (s) does constitute the mirrored counterpart of the characteristic polynomial of the matrix characterizing LM in z-coordinates, once u f is replaced by ρ. After some tedious computations, the perturbed triangular system:

ż1 = z 2 ż2 = z 3 - ω s H e p ż3 = -ω s β 0 H -1 µ v z 1 -(β 1 DH -1 µ v + ω s H -1 µ ω )z 2 -(λ ω + β 1 µ v )z 3 + ω s H -1 (λ ω + β 1 µ v )e p +ω s H -1 µ v η z (z 1 , z 2 , z 3 , e p ) (18)
is derived, where η z (z 1 , z 2 , z 3 , e p ) η V ( δ, ω, Pm , e p ). In accordance with the companion form structure of the matrix characterizing the above error subsystem, it is possible to choose µ v , µ ω , k pω (affecting λ ω ) in order to robustly place all its eigenvalues -namely, the roots of the polynomial Π A (s) -in the left half complex plane C -(independently of X d , T d0 , K c ). Meanwhile, notice how µ ω is to be necessarily positive when D is zero. On the other hand, the e p -dynamics can be derived as:

ėp = -θ 1 P e + θ 2 ω sin 2 (δ) + P e ωcotg(δ) + θ 3 sin(δ)u f -k pω z 3 -ω s H -1 e p + µ v β 0 z 1 -µ v β 1 e p +µ v β 1 Hω -1 s z 3 + µ v β 1 Hω -1 s λ ω z 2 -µ v β 1 k pω z 2 -µ v η z (z 1 , z 2 , z 3 , e p ) + µ ω z 2 , (19)
with P e again satisfying P e = e p -Pm + k pω ω + P m . Finally get:

ėp = -θ 1 e p -θ 1 P m + θ 1 Hω -1 s (z 3 + λ ω z 2 ) -θ 1 k pω z 2 +θ 2 z 2 sin 2 (δ) + P e z 2 cotg(δ) + θ 3 sin(δ)u f -µ v β 0 z 1 -µ v β 1 Dω -1 s z 2 -µ ω z 2 -k pω + µ v β 1 Hω -1 s z 3 + k pω ω s H -1 + µ v β 1 e p +µ v η z (z 1 , z 2 , z 3 , e p ), (20) whereas, from (15)
, write

u f sin(δ) = -k p e p -k r Φ(δ, P e )e p + λ ρ + θ 1 P m θ 3 -k r P ess Λ p , with Λ p = sin 2 (δ) 2 -sin 2 (δ s ) 2 + (P e cotg(δ)) 2 - (P m cotg(δ s )) 2 .
On the basis of the Λ p -expression and the P ess -definition, write

-k r P ess Λ p = k r θ 1 P m θ 3 [k p + k r Φ(δ s , P m )] Λ p α p Λ p , (21)
with α p belonging to the open connected set C p = 0, θ 1 P m θ -1 3 Φ(δ s , P m ) -1 . Now, in the RHS of [START_REF] Marino | A global state feedback output regulating control for uncertain systems in strict feedback form[END_REF], let: A z be the matrix multiplying the vector z = [z 1 , z 2 , z 3 ]

T ; B p be the column vector multiplying e p ; B η be the column vector multiplying η z (z 1 , z 2 , z 3 , e p ). Set the control parameters to satisfy the related assumption of Theorem 2 (concerning the roots of Π A (s)) and let P z be the symmetric positive definite matrix satisfying the Lyapunov equation P z A z + A T z P z = -I in terms of the identity matrix I ∈ M(3, R). Furthermore, consider the composite Lyapunov function

W = z T P z z + 1 2 e 2 p ( 22
)
that satisfies, along the trajectories of the closed loop system, the equality:

Ẇ = -z 2 + 2z T P z B p e p + 2z T P z B η η z (z 1 , z 2 , z 3 , e p ) -θ 1 + k p θ 3 -k pω ω s H -1 -µ v β
e p + θ 1 Dω -1 s -µ ω -µ v β 1 Dω -1 s z 2 e p + θ 1 Hω -1 s -k pω -µ v β 1 Hω -1 s z 3 e p . ( 23
)
Recall that a

λ ρ = -λµ v -β 0 z 1 + β 1 k pω z 2 -β 1 Hω -1 s (z 3 + λ ω z 2 )
+β 1 e p + η z (z 1 , z 2 , z 3 , e p ) + λµ ω z 2 .

Equation ( 23) thus leads to the inequality3

Ẇ ≤ -z 2 + 2 z P z B p |e p | +2 z P z B η |η z (z 1 , z 2 , z 3 , e p )| -θ 1 + k p θ 3 -k pω ω s H -1 -µ v β 1 + λµ v β 1 θ
| + µ v β 0 |z 1 ||e p | +µ v β 0 λθ 3 |z 1 ||e p | + θ 1 Dω -1 s -µ ω -µ v β 1 Dω -1 s +µ v β 1 Dλθ 3 ω -1 s + λθ 3 µ ω |z 2 ||e p | + θ 1 Hω -1 s -k pω -µ v β 1 Hω -1 s +µ v β 1 Hω -1 s λθ 3 |z 3 ||e p | +µ v θ 3 λ|η z (z 1 , z 2 , z 3 , e p )||e p |.
By finally completing the squares and by setting sufficiently large k p and k r (under given µ ω , µ v , k pω ), we obtain, in terms of suitable positive reals g p and g pb , the last compact inequality:

Ẇ ≤ -z 2 /4 -g p e 2 p + g pb |η z (z 1 , z 2 , z 3 , e p )|(|e p | + z ), (24)
where η z (z 1 , z 2 , z 3 , e p ) = h V ( δ, Pe ) just collects the higher order terms of Ṽt in [START_REF] Marino | Learning control for nonlinear systems in output feedback form[END_REF]. Such a Ṽt is injected into the error dynamics by the control [START_REF] Mareels | A simple selftuning controller for stably invertible systems[END_REF], through the gain µ v . Exponential convergence of the (z, e p )-errors is thus achieved for sufficiently small initial (z, e p )-errors: z 1 (0) = δ(0), z 2 (0) = ω(0), z 3 (0) = -λ ω ω(0)+ω s H -1 (P m -ρ(0)-P ess ), e p (0) = P e (0) -P m -k pω ω(0) + (P m -ρ(0) -P ess ), that make, for any t and in accordance with (24): i) Ẇ be negative definite; ii) generator states and references belong to B. With this respect, recall that: i) γ i (i = 1, 2, 3, 4) are positive reals that do not increase while k p and k r increase; ii) P ess turns out to decrease as k p and k r increase. On the other hand, achieving exponential convergence of the (z, e p )-errors to zero guarantees exponential convergence of δ(t) = δ(t)-δ s , ω(t) = ω(t) to zero, along with exponential convergence to zero of Pe (t) and P m -ρ(t) -P ess .

Comment C. The design of ( 15) and the corresponding stability analysis have taken advantage from the y-based arguments of Section III, as the Pm -dynamics has shown. On

V. SIMULATION RESULTS IN SIMSCAPE ELECTRICAL

This section illustrates the previous theoretical derivations concerning the output feedback version of the control algorithm [START_REF] Mareels | A simple selftuning controller for stably invertible systems[END_REF], in comparison to the standard PSS-AVR structure (see Section II.C). The quantities X s and δ required by [START_REF] Mareels | A simple selftuning controller for stably invertible systems[END_REF] have been computed via the (δ, X s )-formulas of [START_REF] Damm | Adaptive nonlinear output feedback for transient stabilization and voltage regulation of power generators with unknown parameters[END_REF], here using X d and just the nominal value of the infinite bus voltage V s . The simulations have been carried out with reference to the Matlab/Simulink Simscape Electrical benchmark (including detailed electrical models and several additional elements), in which the duration of the short circuit -placed just after the transformer connecting the test machine to the transmission line -has been set to 90 ms and no SVC has been considered. In particular, simulations concern the Single Machine Infinite Bus real case at https://fr.mathworks.com/help/physmod/sps/examples/svc-andpss-phasor-model.html, represented by a 1000 MW machine, connected to a 5000 MW power plant through a 700 km line. The physical saturation for the generators exciters, which is an important limitation for the control performance, is the same for both the controllers. The parameters for the standard PSS-AVR are set as: gain= 2, wash-out time constant 0.7, lead-lag time constants T num = 60e -3 , T den = 0.5. The parameters of the output feedback control [START_REF] Mareels | A simple selftuning controller for stably invertible systems[END_REF] are chosen as: k pω = 70, λ = 90, k p = 10, k r = 10, µ v = 10, µ ω = 10. Figure 3 shows how the standard PSS-AVR is not able to keep the system stability: δ goes outbounds. On the other hand, when [START_REF] Mareels | A simple selftuning controller for stably invertible systems[END_REF] is used in the test machine -with the infinite bus remaining with the standard PSS-AVR -, systems stability is assured even in the presence of 5% noise5 on the active and reactive electrical powers P e and Q, as shown in Figure 4: all states are kept inside their operation bounds. According to Figure 5, the standard PSS-AVR input saturates while being not able to bring the generator back to synchronism to the main grid, whereas the (15)-control input saturates as well but it is then able to bring the generator back to its equilibrium point. have been first derived in Section III (see Lemma 1, Theorem 1 and Corollaries 1-2). Such new interpretations have led to the nonlinear control [START_REF] Mareels | A simple selftuning controller for stably invertible systems[END_REF], whose linear action coincides with the one provided by (3): simplicity of control structure (just one integrator is involved) and robust tuning procedure of the linear design (see Theorem 2) have been preserved, with, additionally, no use of the mechanical input power. Control [START_REF] Mareels | A simple selftuning controller for stably invertible systems[END_REF] has endowed the robust coordinated PSS-AVR (3) with: i) a nonlinear term that relies on the power angle, allowing for an automatic gain scheduling; ii) nonlinear robustifying terms on the back-stepping-based electrical power tracking error, with the aim of enlarging the machine stability region (as illustrated by numerical analysis). Realistic simulations concerning a benchmark scenario have finally shown that the nonlinear controller preserves the stability when a mslong short circuit occurs, whereas instability arises when the standard PSS-AVR is used.

2 -b 2 =

 22 (a -b)(a + b) for a, b ∈ R; then write P e cotg(δ) -P m cotg(δ s ) = P e cotg(δ) -P m cotg(δ) + P m cotg(δ) -P m cotg(δ s ), as well as Pe = e p + k pω z 2 -Hω -1 s (z 3 + λ ω z 2 ). Then obtain |Λ p | ≤ 3 i=1 γ pi |z i | + γ 4 |e p | over B ⊂ D, with γ i (i = 1, 2,[START_REF] Bourlés | Analysis and design of a robust coordinated AVR/PSS[END_REF][START_REF] Damm | Transient stabilization and voltage regulation of a synchronous generator[END_REF] being positive reals depending on the diameter of B but not increasing while k p and k r increase. On the other hand, write

Fig. 1 .

 1 Fig. 1. Power angle δ, relative speed ω, electrical power Pe, terminal voltage Vt (blue: DFLR; red: proposed nonlinear controller).

Fig. 2 .

 2 Fig. 2. Power angle δ, relative speed ω, electrical power Pe, terminal voltage Vt (blue: DFLR (3); red: nonlinear controller (15)).

Fig. 3 .

 3 Fig. 3. Standard PSS-AVR: (a) power angle δ in degrees; (b) angular speeds ω 1,2 in p.u.; (c) output voltages V t1,2 .

Fig. 4 .

 4 Fig. 4. Output feedback control (15) under noise effects: (a) power angle δ in degrees; (b) angular speeds ω 1,2 in p.u.; (c) output voltages V t1,2 in p.u.

Fig. 5 .

 5 Fig. 5. Standard PSS-AVR (left); output feedback control (15) under noise effects (right): control input in p.u.

  1 e 2 p +α p Λ p θ 3 e p + λ ρθ 3 e p -k r θ 3 Φ(δ, P e )e 2 p +θ 2 sin 2 (δ)z 2 e p + P e cotg(δ)z 2 e p +µ v η z (z 1 , z 2 , z 3 , e p )e p -µ v β 0 z 1

  3 e 2 p +α p |Λ p |θ 3 |e p | -k r θ 3 Φ(δ, P e )e 2 p +θ 2 sin 2 (δ)|z 2 ||e p | + |P e cotg(δ)||z 2 ||e p | +µ v |η z (z 1 , z 2 , z 3 , e p )||e p

Uncertainties about the infinite bus voltage value are usually present in the practice. In the realistic simulations of Section V, just the nominal value of the infinite bus voltage will be used.

Authorized licensed use limited to: IFSTTAR. Downloaded on March 05,2021 at 09:44:04 UTC from IEEE Xplore. Restrictions apply.

Notice that, when λθ 3 = 1, several coefficients disappear from the nondefinite sign terms. This explains the feedforward role of the λ-term in[START_REF] Mareels | A simple selftuning controller for stably invertible systems[END_REF].

The use of the proposed technique can be extended to similar scenarios concerning nonlinear systems with extended matching uncertainties.Authorized licensed use limited to: IFSTTAR. Downloaded on March 05,2021 at 09:44:04 UTC from IEEE Xplore. Restrictions apply.

Usually (in similar scenarios) all current structures work pretty well under noise. Authorized licensed use limited to: IFSTTAR. Downloaded on March 05,2021 at 09:44:04 UTC from IEEE Xplore. Restrictions apply.

the other hand, a back-stepping action that is reminiscent of the typical two-time-scale arguments has been adopted 4 . It is confined to the natural definition of the electrical power reference, so that all the uncertainties in the corresponding electrical power tracking error dynamics are compensated by the integrator placed into the upper mechanical subsystem. Nonlinear robustifying actions can be just included at last, through the use of a composite Lyapunov function. In contrast to the fully linearized analyses which Theorem 1 and Corollary 2 rely on, genuinely nonlinear stability tools (applied to the nonlinear model [START_REF] Anderson | Power System Control and Stability[END_REF] with just local validity) allow us to neglect, in the Lyapunov analysis of Theorem 2, just the nonlinearities concerning the terminal voltage equation, whereas the well-known nonlinearities affecting the electrical power dynamics are explicitly faced by the nonlinear features of the controller (aiming at enlarging the stability region) and by one last control parameter to be set at last.

Comment D. According to Comment B, set µ ω = 0, ρ(0) = P m , replace δ appearing at the denominator by δ s and neglect higher order robustifying actions in [START_REF] Mareels | A simple selftuning controller for stably invertible systems[END_REF]. Then, (15) coincides with DFLR (3) [under equivalent gains], which admits an attraction domain R l (whose estimate may be provided by linearization arguments), whereas a possibly different attraction domain R nl exists (whose estimate is given by the proof of Theorem 2) for the corresponding nonlinear version [START_REF] Mareels | A simple selftuning controller for stably invertible systems[END_REF]. The effectiveness of [START_REF] Mareels | A simple selftuning controller for stably invertible systems[END_REF] in including points in the state space of (1) -that do not belong to R l -into its attraction domain R nl can be inspected through a numerical analysis. Consider the synchronous generator in [START_REF] Verrelli | Output feedback transient stabilization and voltage regulation of a synchronous generator[END_REF] (no input saturations are considered), whose physical parameter values are: ω s = 314.159 rad/s, D = 5 p.u., H = 8 s, T d0 = 6.9 s, K c = 1, X d = 1.863 p.u., X d = 0.257 p.u., X T = 0.127 p.u., X L = 0.4853 p.u.. The nominal operating condition δ s =1.256 rad, P m =0.9 p.u., ω=0 rad/s guarantees V t = V tr = 1 p.u., with V s = 1 p.u.. The common control parameters for both the controllers ( 15) and (3) [under equivalent gains] are (in SI units): k pω = 1.5, λ = 3, k p = 1.5, k r = 9, µ v = 18 (µ ω = 0), while the controller initial conditions are compatible with the operating point (with the feedforward contribution u f s being included in the controllers). Figure 1 concerning the numerical integration of the equations involved in the considered scenario -shows that there exists a set of initial conditions, namely δ(0) = 0.3 rad, ω(0) = -3 rad/s, P e (0) = 2.2 p.u., that belongs to the stability region R nl guaranteed by [START_REF] Mareels | A simple selftuning controller for stably invertible systems[END_REF] but not on R l .

Comment E. The nonlinear robust coordinated PSS-AVR (15) can be seen as the evolution of the nonlinear decentralized robust adaptive controller in [START_REF] Verrelli | Robust transient stabilization problem for a synchronous generator in a power network[END_REF], which was exclusively designed for transient stabilization, under a rather weak set of assumptions on the remote network dynamics. The inclusion in [START_REF] Mareels | A simple selftuning controller for stably invertible systems[END_REF] of the terminal voltage regulation error to drive the mechanical input power estimate (acting as the integral term of (3)) allows [START_REF] Mareels | A simple selftuning controller for stably invertible systems[END_REF] to achieve voltage regulation even at permanently perturbed operating points (see the subsequent Figure 2 concerning a 70%-decrease of P m and a double


Yu Guo
email: yu.guo@etudiant.univ-reims.fr

Su Ruan
email: su.ruan@univ-reims.fr

Jérôme Landré

Jean-Marc Constans
email: constans-jm@chu-caen.fr

Magnetic Resonance Spectrum Separation Using Sparse Representations and Wavelet Filters

Magnetic Resonance spectroscopy (MRS) provides a "frequency-signal intensity" spectrum of multiple peaks that reflect the biochemical composition of a localized region in the body. The peak intensity or the area under each peak is proportional to the concentration of that assigned metabolite. Accurate quantification of in vivo MRS (measuring peak intensity or area) is very important to diagnose certain metabolic disorders. However, strongly overlapping metabolite peaks, poor knowledge about background component (the baseline), and low signalto-noise ratio (SNR) make the task difficult. In this paper, a novel spectrum separation method using sparse representations and wavelet filters is proposed to separate baseline and spectra of different metabolites and finally achieves an accurate MRS quantification. With the proposed method, the accuracy and the robustness of MRS quantification are improved, from simulation data, compared with a commonly used frequency-domain MRS quantification method. The quantification on tumor metabolism with in vivo brain MR spectra is also demonstrated.

I. INTRODUCTION

Over the years MRS (Magnetic resonance spectroscopy) has clearly proven to be an important non-invasive tool in biomedical research to study humans and animals [START_REF] Hja In't Zandt | Common processing of in vivo MR spectra[END_REF]. The observed MRS spectrum is a mixture of several spectra respectively corresponding to different metabolites in a localized region under study. Each biochemical or metabolite is identified by its unique position or chemical shift along the frequency axis of a spectrum [START_REF] Stanley | Magnetic Resonance Spectroscopy in Schizophrenia: Methodological Issues and Findings-Part I[END_REF]. The biochemical quantities (concentration, PH) of a certain metabolite can be obtained by accurately quantifying the corresponding MR spectrum (computing the peak intensity or the area under the spectrum peak). Accurate quantification of in vivo MRS is very important to diagnose certain metabolic disorders, especially those affecting the brain. However, when a large background component (the baseline) is present along with a low signal to noise ratio, accurate quantification of MR spectra is a challenging problem [START_REF] Bartha | Factors affecting the quantification of short echo in-vivo 1H MR spectra: prior knowledge, peak elimination, and filtering[END_REF]. In this paper, a method using sparse representations and wavelet filters is proposed to separate baseline and spectra of different metabolites for achieving more accurate MRS quantification.

In recent years sparse representation is attracting more and more interest in signal processing domain. In the study of sparse representation, in literature [START_REF] Mallat | Matching pursuits with time-frequency dictionaries[END_REF]- [START_REF] Rao | Subset selection in noise based on diversity measure minimization[END_REF] lots of pursuits algorithms are proposed for finding the sparsest representation of a signal in a given dictionary. Taking into account the a priori knowledge about the model of MR spectra and the peak frequencies of different metabolites, the representation of a mixed MR spectrum in a dictionary constructed with a set of Lorentizian and Gaussian basis functions will be sparse and disjoint. Therefore, a non-negative pursuit algorithm which is based on the regulative FOCUSS algorithm in [START_REF] Rao | Subset selection in noise based on diversity measure minimization[END_REF] is used to estimate the sparse representation and separate spectra of different metabolites.

Besides, for eliminating the influence of the baseline to the estimation of sparse representation coefficients, a wavelet filter is used to remove the smooth component of an observed spectrum. Because of the smoothness and broadness of baseline, baseline will be eliminated with the wavelet filters. However, the overlapping component of metabolite spectra with the baseline will also be removed with the filter. By processing each basis function in the given dictionary with the same filter, the estimation of representation coefficients will not be influenced by the removed component of metabolite spectra.

The proposed method is evaluated with simulation experiments and compared with a commonly used frequencydomain MRS quantification method in [START_REF] Gillies | Quantification of MRS data in the frequency domain using a wavelet filter, an approximated Voigt lineshape model and prior knowledge[END_REF]. The results of these experiments show the effective performance of the proposed method. The quantifications of in vivo brain MR spectra are confirmed by the experts and demonstrated in this paper.

II. THEORY AND METHOD A. Signal model

Generally, the mathematical model of a mixed metabolite spectrum S can be expressed in the following form:

S = K k=1 s i = K k=1 [L k (f ) + G k (f )] (1)
with

L k (f ) = a Lk 1 + ((f -f Lk )/d 2
Lk ) and

G k (f ) = a Gk e -((f -f Gk )/d Gk ) 2
where L k and G k denote Lorentzian and Gaussian lineshape respectively, f is the frequency of each data point, K stands for the number of lines used to build up the spectrum and the spectrum of single metabolite is corresponding to one or several lines [START_REF] Mierisova | MR spectroscopy quantitation: a review of frequency domain methods[END_REF]. p k = [a Lk , d Lk , f Lk , a Gk , d Gk , f Gk ] is a model parameter vector, where a is the intensity, d the linewidth and f the central frequency. The spectra of different metabolites have different central frequencies which are known and can be used as an a priori knowledge in a MRS quantification method. When a possible baseline contribution and a noise which is often assumed as Gaussian distributed are considered, an observed MR spectrum can be modelized as

X = S + B + e (2)
In most of MRS quantification approaches, especially the data obtained by quantum mechanics or from measured in vitro metabolite solutions are unavailable such as the method in [START_REF] Vanhamme | Improved method for accurate and efficient quantification of MRS data with use of prior knowledge[END_REF], quantification can be achieved by using a nonlinear leastsquares algorithm to estimate the nonlinear parameter vector p. To deal with the baseline problem, baseline is generally approximated by a mathematical function (a sum of splines, wavelets, sinusoids or polynomials with adjustable parameters) then included in the parametric nonlinear least-squares fit or separately accommodated in a pre-processing step [START_REF] Ratiney | Time-domain semi-parametric estimation based on a metabolite basis set[END_REF].

For these MRS quantification methods based on nonlinear optimization, some nonlinear methods have the risk of converging to a local minimum while the other have the disadvantage of poor computational efficiency, and their performances deteriorate with the increase of the number of model parameters. For a better application, the quantification accuracy needs to be improved, especially in the condition that large baseline components overlap with metabolite spectra. For avoiding the disadvantage caused by nonlinear optimization and improving the quantification accuracy, a linear non-parameter method proposed in this paper.

B. Sparse representation

A sparse representation model can be described as

x = Dw (3)
where x N ×1 is the signal vector which can be sparsely represented as a linear combination of columns (often called "basis vectors") of a dictionary matrix

D N ×M = [d 1 , d 2 , . . . , d M ],
and when M > N , it is a overcomplete dictionary. w M ×1 is the representation coefficient vector, minimal elements of which are nonzero.

The study of sparse representation has concentrated mainly on the design of dictionaries and pursuit algorithms. Designing dictionaries to better fit the sparse representation model can be done by either selecting one from a prespecified set of linear transforms or adapting the dictionary to a set of training signals [START_REF] Aharon | An Algorithm for Designing Overcomplete Dictionaries for Sparse Representation[END_REF]. When the dictionary is known and fixed, pursuit algorithms can be used for the extraction of the sparsest representation. The main pursuit algorithms in literature are matching pursuit [START_REF] Mallat | Matching pursuits with time-frequency dictionaries[END_REF] and the orthogonal matching pursuit (OMP) algorithms [START_REF] Pati | Orthogonal matching pursuit: Recursive function approximation with applications to wavelet decomposition[END_REF] which are greedy algorithms that select basis functions sequentially by involving the computation of inner products between the signal and dictionary columns, the basis pursuit algorithm which suggests a convexification of the problems by replacing l 0 -norm with l 1 -norm [START_REF] Chen | Atomic decomposition by basis pursuit[END_REF], the focal underdetermined system solver (FOCUSS) which using the l p -norm with p ≤ 1 as a replacement for the l 1 -norm [START_REF] Gorodnitsky | Sparse signal reconstruction from limited data using FOCUSS: A re-weighted norm minimization algorithm[END_REF][START_REF] Rao | Subset selection in noise based on diversity measure minimization[END_REF].

In this paper sparse representation is used for separating MR spectra. Based on the character and the available a priori knowledge of the particular signal, the dictionary is designed and a non-negative pursuit algorithm is used to finding the sparsest representation.

C. Dictionary construction

According to the model of MR spectra in (1), a mixed metabolite spectrum can be represented as the linear combination of several Gaussian and Lorentizian functions. Hence, an overcomplete dictionary D, which contains enough basis vectors designed as the sample vectors of normalized Gaussian and Lorentizian functions, can will fit the sparse representation model in [START_REF] Bartha | Factors affecting the quantification of short echo in-vivo 1H MR spectra: prior knowledge, peak elimination, and filtering[END_REF]. Then, a mixed metabolite spectrum vector S can be represented as

S = Dw = M m=1 d m w m
where d m is one of the basis vectors which has the same length as S, M is the number of basis vectors in D. For a detailed illumination of the construction of the dictionary, all the basis vectors in D are divided into K groups {D 1 , D 2 , . . . , D K }, where K is the number of the spectrum peaks under study. The basis vectors d ij (j = 1 . . . L) in group D i (i = 1 . . . K) have the same central frequencies with the known central frequency of one of the metabolite peaks under study and different widths d ij = d 0j +j∆d. Because the spectrum peaks of different metabolites have different central frequencies, the representations of spectra of different metabolites in the constructed dictionary are disjoint. If we can find w the sparsest solution of S = Dw, we can decompose the spectrum S into several components S i = D i w i (i = 1 . . . K) which relate to spectra of different metabolites.

However, because of the presence of baseline B in an observed spectrum, the sparest solution w of X = Dw will certainly depart from the expected solution of S = Dw. The smoothness of baseline is usually used to deal with baseline problem and a wavelet filter can be used to eliminate baseline [START_REF] Ratiney | Time-domain semi-parametric estimation based on a metabolite basis set[END_REF]. However, because of the overlapping between baseline and the metabolite spectra, when a wavelet filter is used to process observed spectra, it is difficult to only remove baseline and reserve the whole mixed metabolite spectrum.

Here, for accurate estimation of the sparest solution w, each basis vector of the constructed dictionary is also processed with the same wavelet filter, then a basis vector d in D can be denoted as d = d h + d l , where d l is the removed part of d with the wavelet filter, and d h the reserved parts. Based on the wavelet filter processing results of all the basis vectors in D, the dictionary D can be divided into D h and D l , which are corresponding to d h and d l respectively. The reserved part S h of a mixed metabolite spectrum S with the wavelet filter can be represented as S h = D h w. When the wavelet filter removes the whole baseline, there is S h = x h , where x h the reserved part of an observed spectrum x. Furthermore, if each column of D h is not zeros or too small value vector, the sparest solution w of x h = D h w is equal to the sparse representation of a mixed metabolite spectrum S in the constructed dictionary D.

Therefore, MR spectra separation can be achieved by finding the sparsest solution w of x h = D h w.

D. Non-negative FOCUSS pursuit algorithm

Theoretically, the sparsest representation of y N ×1 with respect to an overcomplete dictionary D N ×M is the following optimization problem: where . 0 is the l 0 norm, for counting the nonzero entries of a vector; w M ×1 is the coefficient vector of basic functions. However, it is a NP-hard problem. Thus, approximate solutions are considered instead, the FOCUSS minimizes l (p≤1) -norm in place of l 0 -norm to obtain sparse solution. Then, the sparse representation becomes the solution of

min w M i=1 sgn(w(i)) | w(i) | p subject to y = Dw (6)
When the noise exists, an exact minimum norm solution can not be sought. Instead, a solution that minimizes l (p≤1) -norm and approximately satisfies the set of constraints is found. The solution is

w = arg min w J(w) (7)
where

J(w) = [ Dw -y 2 + γE (p) (w)]
For the problem in this paper, taking into account of the non-negativity of the representation coefficients, a nonnegative constraint is added. So the representation coefficient vector w can be estimated as

w = arg min w J(w) (8)
where

J(w) = [ Dw -y 2 + γE (p) (w)] and ∀i : w i ≥ 0
At each iteration step, we set the negative values of the solution to zero for guaranteeing the nonnegative character of coefficient vector w. Based on the basic iterative form of the regularized FOCUSS algorithm, the iterative form of the nonnegative pursuit algorithm becomes:

W k+1 = diag(w 1-(p/2) ki ) (9a) w k+1 = W k+1 D T k+1 (D k+1 D T k+1 + λI) -1 y (9b) w k+1 (i) = 0 if w k+1 (i) < 0 w k+1 (i) if w k+1 (i) ≥ 0 (9c)
The parameter λ controls the tradeoff between quality of fit y -Dw and the degree of sparsity, and the value of λ should increase with the level of noise. The iterative algorithm in ( 9) is used to estimate the sparsest solution w of x h = D h w. With the representation coefficient vector w, spectrum of different metabolites can be separated by computing s i = D i x i (i = 1 . . . K). Finally, MRS quantification can be achieved easily benefiting from the separation results.

III. EVALUATION AND RESULT

A. Evaluation with simulated spectra

The proposed approach is tested with simulated 1 H human brain MRS signals. A simulated MR spectrum consisting of 512 data points is composed of eleven single metabolite spectra modelized as Gaussian functions, a baseline mixing of nine models and a certain level of Gaussian noise. The parameters of simulated metabolite spectra and baseline are shown in Table I. For each noise level, a set of 100 simulated MR spectra is created in order to give reliable values for the uncertainties of the estimation results and check the robustness of the proposed method. The value of the SNR is defined as the ratio of the amplitude of the NAA resonance to the standard deviation of noise e.

As an a priori knowledge, the basis vectors in the constructed dictionary supposed to represent spectra of main metabolites under study have the known central frequencies as shown in Table I and the unknown linewidths taken from the range [d 1 , d 2 ] = [0.5, 10] with the sample step ∆ d = 0.5. In this paper, the Coiflet discrete wavelet decomposition of the observed MR spectra using COIF5 at level 5, is computed and the detail coefficients of the wavelet decomposition are reserved to reconstruct x h , and the same wavelet decomposition

B. Results

Figure 1 displays a simulated MR spectrum with SNR=17dB as in vivo conditions, and the corresponding separation results. The comparison of separation results and real spectra in Fig. 1 shows that baseline is well separated from metabolite spectra of interest.

Figure 2 plots the RRMSE (Relative Root Mean Square Error) of estimated metabolite peak areas in different noise levels (SNR =12dB, 14dB, 16dB, 18dB, 20dB). RRMSE is defined as the ratio of RMSE (Root Mean Square Error) of estimated results to the real value. Because of the different peak intensities, a noise influences the quantification of these metabolites in different degrees. As demonstrated in Fig. 2, for the metabolites with relatively big peaks, such as NAA, Cr, the estimation results are less influenced with the increase of the noise, while for the weak metabolite peaks, especially the metabolites Lac and Cho, the noise have relatively big influence to the estimation result. Like most of the other MRS quantification methods, the estimation accuracy decreases with the increase of noise level.

Using simulated MR spectra with SNR=17dB, the proposed method is compared with a commonly used frequency-domain MRS quantification method in which Levenberg-Marquardt algorithm is used to estimate the nonlinear model parameters of metabolite spectra and a wavelet filter is used to remove the baseline component in an iterative subtraction manner [START_REF] Gillies | Quantification of MRS data in the frequency domain using a wavelet filter, an approximated Voigt lineshape model and prior knowledge[END_REF]. In simulation experiments, the same a priori knowledge is used in both methods. For the method in [START_REF] Gillies | Quantification of MRS data in the frequency domain using a wavelet filter, an approximated Voigt lineshape model and prior knowledge[END_REF], the choice of the initial parameters can deeply influence the estimation result, while in the proposed method, a small random positive value vector is used for the initialization of the representation coefficients and different initializations have no influence to estimation results. The RRMSE of metabolite peak areas and the error bars (the standard deviation of relative estimation errors) of the both methods are illustrated in Fig. 3. Obviously, the relatively small estimation errors and error bars of the proposed method show the superiority of our method.

C. In vivo MRS data

The method was also tested on in vivo 1 H human brain MR spectra. Two human brain MR spectra from a brain tumor patient obtained with PRESS and an echo-time of 35 ms were quantified by our method.

IV. CONCLUSIONS

The new MRS quantification method proposed in this paper combines a wavelet filter and a linear non-parametric algorithm based on sparse representation to remove baseline and separate spectra of different metabolites and finally achieves an accurate MRS quantification. The method can separate the overlapping components of baseline and metabolite spectra of interest and well deals with the challenging baseline problem.

In addition, it effectively uses the a priori knowledge and achieves satisfying MR spectrum separation. Compared with commonly used nonlinear parametric optimization algorithms, the proposed method has a better computational efficiency, avoid the choice of initial parameters and are more robust. In future work, the method will be tested on large data of in vivo brain spectra.

  to y -Dw 2 ≤ (5)

Fig. 1 .

 1 Fig. 1. Separation results of a simulated MR spectrum (SNR=17dB) using the proposed method (a) raw spectrum; (b) the estimated and the true baseline (dashed line); (c) the estimated and the true spectrum (dashed line).

Fig. 2 .

 2 Fig. 2. The RRMSE of estimated metabolite peak areas in different noise levels with the proposed method.

Fig. 3 .

 3 Fig. 3. Comparison results. RRMSE of estimated metabolite peak areas and the error bars of the proposed method and the nonlinear one (SNR=17dB).

Fig. 4 (

 4 a) and Fig. 4(c) are the MR spectra of tumor tissue and normal tissue respectively. Fig. 4 (b) and Fig. 4 (d) show the separation results of the baseline and metabolite spectra of interest. According to the spectrum separation results, the peak areas of the main metabolites in each MR spectra were computed. Comparing the different MRS quantification results of normal tissue and tumor tissue, doctors can analyse the kind of tumor more easily.

Fig. 4 .

 4 Fig. 4. The separation results of in vivo 1H human brain MR spectra. (a) spectrum from a tumor tissue; (b) separation results of the MR spectrum in (a); (c) spectrum from a normal tissue; (d) separation results of the MR spectrum in (c) (BL=baseline, MSI=spectra of interest).

TABLE I PARAMETERS

 I OF SIMULATED METABOLITE SPECTRA AND BASELINE.

	Parameters of metabolic spectra		Parameters of baseline
	Metab.	f k (ppm)	d k	a k	k	f k (ppm)	d k	a k
	Cr	3.91	8.0	8.0	1 G	2.02	85	2.33
	Glu/Gln	3.74	4.0	4.0	2 G	2.35	100 0.33
	mI	3.56	4.0	5.5	3 G	2.50	70	0.67
	Tau	3.42	3.0	3.2	4 G	3.00	70	0.67
	Cho	3.22	6.0	2.0	5 G	3.29	100	1.0
	Cr/PCr	3.03	6.0	4.0	6 G	3.50	80	2.00
	GABA	2.37	4.0	1.5	7 G	4.00	70	0.67
	Glu/ Gln	2.12	6.0	9.5	8 L	0.90	300 1.33
	NAA	2.02	5.0 13.0	9 L	1.30	300 1.67
	Lac	1.33	5.0	1.5
	Lac	1.26	3.0	1.5