Bernard Xerri
email: bernard.xerri@univ-tln.fr

Bruno Borloz
email: bruno.borloz@univ-tln.fr

The Multicriteria Constrained Stochastic Matched Filter For Underwater Bioacoustic Signals

Keywords: detection, bioacoustics, Signal to Noise Ratio (SNR), stochastic matched filter, multi-noise, correlation matrix

The aim of this paper is the detection of a bioacoustic signal embedded in several noises such as sea noise and other bioacoustic signals (dolphins, sperm whales). All the signals are real world signals. Only second order statistics are use through the estimated correlation matrices of the signals. This paper proposes an extension of the Constrained Stochastic Matched Filter (CSMF) based on the optimization of the Signal to Noise Ratio after linear filtering. The approach proposed is a multicriteria one, merging three different versions of the CSMF, and is named Multicriteria CSMF (MCSMF).

The objective is that the results obtained are better than the other methods, or at least equal to the best among the three.

The results are provided on ROC curves and the method is compared to the classical method Stochastic Matched Filter (SMF).

I. INTRODUCTION

Detection of bioacoustic signals is an important area of research to, for instance, estimate animal density or abundance from acoustic data.

In this paper we consider dolphin and sperm whales signals embedded in sea noise : these signals are real world ones found in the SABIOD/Bombyx Database [START_REF]SABIOD Scaled Acoustic[END_REF]. Dolphins can produce signals such as whistles or clicks. We are interested only in clicks (more difficult to detect). Our objective is to detect dolphins signals, without confusion with the others, with a linear filter which is an extension of the Stochastic Matched Filter (SMF), named Multicriteria Constrained Stochastic Matched Filter (MCSMF).

A test signal is created with realizations of dolphins and sperm whales clicks with a normalized energy, and sea noise the power of which can be modified to impose the signal to noise ratio (SNRsea). The sea noise is always there and the clicks are added up.

The method we propose is compared to two other classical method: the first one is the SMF for which the noise is only the sea noise. The second one is the SMF for which the noise is the mean of the undesired signals (sea noise and sperm whale signals).

All the methods tested only use correlation matrices of the different signals.

II. THE MULTICRITERIA CONSTRAINED STOCHASTIC MATCHED FILTER

We consider that the × -dimension correlation matrices of each class of signals are known or estimated from data learning.

Here we try to detect dolphin clicks (correlation matrix ) when other signals are unwanted ( , are the correlation matrices for sea noise and sperm whales).

A. The SMF

If the signal to detect were determinist, the optimal filter to use is the Matched Filter (MF) which consists in projecting the measured data onto a one dimensional ( = 1 ) subspace spanned by ( is the noise correlation matrix).

When the signal is stochastic, the question is: "is still equal to 1 for optimal detection ?".

Classically the Stochastic Matched Filter (SMF) is used when there is a useful signal, of correlation matrix , to detect and only one disturbing noise of correlation matrix .

It consists in calculating the eigenelements of , projecting the measured data onto each eigenvectors and summing ≥ 1 powers of these terms. [START_REF] Bouffaut | Passive Stochastic Matched Filter for Antarctic blue whale call detection[END_REF][START_REF] Bouffaut | Antarctic Blue Whale calls detection based on an improved version of the stochastic matched filter[END_REF][START_REF] Courmontagne | The Adaptive Stochastic Matched Filter for SAS Images denoising[END_REF].

The optimal value of is chosen according to the different ROC curves. Nevertheless, the obtained dimensional subspace has no reason to be the one where the SNR is maximal.

In a multi noise context, people generally consider a mean noise of correlation matrix = + (which are the correlation matrices of the two undesired noises) and then use the SMF to solve this problem. It will be called SMF Σ in the following.

B. The CSMF

Another approach consists in searching a -dimension projection subspace ( < ) spanned by a × orthonormal matrix , in which the Signal to Noise Ratio (SNR) given by = is maximal.

Again, the optimal value of is chosen according to the different ROC curves, and has no reason to be the same than the SMF's one. This filter is the Constrained Stochastic Matched Filter (CSMF) [START_REF] Borloz | Subspace SNR Maximization: The Constrained Stochastic Matched Filter[END_REF].

This has been extended in a multi noise context [START_REF] Borloz | Estimation, détection, classification par maximisation du rapport signal-à-bruit : le filtre adapté stochastique sous contrainte[END_REF]. Then we can define for each noise a SNR ρ =

. The objective is

to maximize a criterium like ! = "∏ $ % $&' or ' % ∑ $ % $&' or % ∑ ) * + , +-)
(geometric, arithmetic or harmonic mean).

All the correlation matrices are full-rank and normalized (trace equal to 1).

Those three criteria lead to a similar equation in which appear and a weighted mean of the . , ∑ α B % &'

(the αi depend of the chosen criterium). These equations are solved with the same algorithm (CSMF) which leads to and the $ (hence !).

Denoting 1' , 1% and 13 ∈ 5 6×1 the × -matrices of the vectors spanning the three optimal subspaces (geometric, arithmetic and harmonic), we calculate the power of the measured signal x after projection onto 1$ , 7 = 1,2,3 :

tr X = > xx > X = = x > X = X = > x
which is used to perform detection (with a threshold).

Note that the harmonic criterium leads to the CSMF solution using = + .

C. The MCSMF

The experiments show that no criterium is always better than the others. Hence, a new idea consists in merging the three optimal subspaces after a power normalization. This is the Multicriteria CSMF.

We define

1@ = A 1' "B ' 1% "B % 13 "B 3 C ∈ 5 6×(31)
with B $ = tr 1$ > 1$ , 7 = 1,2,3. Then we calculate the power of the measured signal F after projection onto 1@ : tr 1@ > FF G 1@ = F G 1@ 1@ > F which is used to perform detection (with a threshold).

The hope is that the results obtained are better than the other methods, or at least equal to the best among the three.

III. EXPERIMENTS

From the SABIOD database, we have extracted a few realizations of bioacoustics signal (dolphins and sperm whales). It is a hard and time-consuming task because the database is poorly informed: hence we have only 15 realizations of clicks for each category.

However, the sea noise is easy to extract and we have a long record available.

A. The bioacoustic signals

Fig. 1 shows one out of the 15 realizations of dolphin clicks: its length is 150 samples (fs=44kHz).

B. The test signal

A test signal is created by addition of different realizations of the bioacoustic signals in sea noise which is present everywhere.

The test signal contains 40 realizations of dolphins clicks and 40 realizations of sperm whales clicks. Necessarily there are repetitions of some realizations, but however the sea noise is each time different and so the situation is never the same.

We don't have the choice, having a few number of realizations and must use some of them several times: this is not a satisfying situation but it is necessary to estimate acceptable ROC curves.

C. Correlation matrices

Considering H realizations . (of length ) of a signal, the correlation matrix I is an × -matrix of rank H obtained by:

I = J . . > K $&'
Unfortunately, we have = 150 but only 15 realizations available: in order not to take all the realizations for the estimation of the correlation matrix and keep some for the test signal, we choose H = 10.

Using the SMF or CSMF theory, usually, correlation matrices are full-rank. Here, this is not the case and this is a serious difficulty: it is necessary to define a signal dolphin subspace which is of dimension H = 10.

D. Signal Processing and Results

The usual processing consists in a mobile linear filtering of a piece of length n of the test signal F: > F where is a × matrix spanning the optimal p-dimension subspace.

> F is a -dimension vector and we calculate its square norm to create a signal to be thresholded and make a decision, called decision signal.

Given the above, a preprocessing stage is necessary: it consists in projecting F onto the signal dolphin subspace: the result is F′. In this subspace, we work with ', O and O of dimension H × H and no more , and of dimension × .

Then we have ′ of dimension H × , and the decision signal is the square norm of ′ > F′.

The correlation matrices ′, O and O are full-rank but not normalized. Hence, the different criteria CSMF (arithmetic, geometric and harmonic) need the $ to be weighted according to the traces of the matrices. Fig. 5 shows an example of the obtained decision signal for the multicriteria filter. Areas where dolphin signal is present are drawn with a dotted square function. . The ROC curves are obtained with a varying threshold (from 0 to the max value). They are given on Fig. 6 and7.

On Fig. 6 are compared the results obtained with the three different CSMF presented above and the Multicriteria method (MCSMF).

The MCSMF is here slightly better than the others (equivalent to the geometric criterium).

On Fig. 7 are compared the results obtained with MCSMF and the SMF methods: the SMF where = is the sea noise correlation matrix (SMF Sea), and the SMF where = +

(SMF Σ).

The MCSMF is here better than the others. Note that for the detection of sperm whales, the MCSMF gives equivalent results than SMF Σ which is usually used.

The calculation of P Q is performed as follows: we count the number of samples above the threshold in the detection areas and divide it by the total number of samples in the detection areas. Such an approach is very penalizing and underestimates the detection probability.

Here below is given an example of ROC curves for different values of p for the geometric CSMF. We can observe that the best result is for = 5. We propose an extension of the Constrained Stochastic Matched Filter (CSMF) based on the optimization of the Signal to Noise Ratio after linear filtering. The approach proposed is a multicriteria one, merging three different versions of the CSMF (arithmetic, geometric and harmonic criteria), and is named Multicriteria CSMF (MCSMF).

The MCSMF has been tested on real world bioacoustic signals and compared to the three versions of the CSMF and to the classical method Stochastic Matched Filter (SMF).

The results are provided on ROC curves. For dolphin detection, they are better for the MCSMF than for the other methods, or at least equal to the best among the three CSMF.

Different powers of the sea noise have been tested (SNR of 0 dB and 6 dB) : obviously, the higher is this power, more difficult it is to discriminate the results of SMF and MCSMF.

For sperm whales detection, all the results are equivalents.

Fig. 1 .

 1 Fig. 1. Dolphin click

Fig. 2

 2 Fig.2shows one out of the 15 realizations of sperm whales clicks: its length is 150 samples (fs=44kHz).

Fig. 2 .

 2 Fig. 2. Sperm Whale click

Fig. 3

 3 Fig. 3 shows the test signal for a SNR (vs. sea noise) of 6 dB: dolphins clicks are present only in the tagged areas.

Fig. 3 .Fig. 4 .

 34 Fig. 3. Test Signal (SNRsea = 6 dB)

Fig. 5 .

 5 Fig. 5. Decision Signal (SNRsea = 6 dB)

Fig. 6 .Fig. 7 .

 67 Fig. 6. ROC curve for CSMF methods (SNRsea = 6 dB)

Fig. 8 .

 8 Fig. 8. ROC curve for MCSMF and SMF methods (SNRsea = 0 dB)

Fig. 9 .

 9 Fig. 9. ROC curves for geometric CSMF (SNRsea = 6 dB)


N Lorrière

G Chabriel

J Barrère

M Pasquinelli

G Pic

N Vannieuwenhuyse

L Escoubas

J.-J Simon

LiFi Reception from Organic Photovoltaic Modules Subject to Additional DC Illuminations and Shading Effects

In this paper, we study the performance of organic photovoltaic (OPV) modules as LiFi receivers in two specific configurations. The PV-based LiFi receiver is first exposed to an additional homogeneous light source with different intensity levels and then the influence of partial lighting is studied (shading effect). In both cases, we compare the sensibility and the cutoff frequency of LiFi transmission when the solar cell is operated either in short-circuit mode (i.e. when it is loaded with an active transimpedance amplifier) or in open circuit mode (i.e. when it is terminated with an high impedance passive load). While the OPV module performance decreases in open-circuit mode as a function of the DC illumination level, we observe an improvement of the cutoff frequency in short circuit mode. This result seems very promising for outdoor LiFi transmissions but also for indoor conditions where natural light can disturb LiFi communications. Theoretical explanations involving physical parameters for energy harvesting (carrier mobility, lumped series resistance) are proposed to justify the observed behaviors. Finally, experimental results of shaded solar cells are provided in the two operating modes (open-circuit and short-circuit). We show that sensibility and bandwidth of OPV modules strongly depend on both the shading configuration and the operating mode.

I. INTRODUCTION

The possibility of using photovoltaic (PV) modules for visible light communications (VLC) is still relatively new [START_REF] Kim | Simultaneous Reception of Visible Light Communication and Optical Energy using a Solar Cell Receiver[END_REF]- [START_REF] Sindhubala | Receiver Intend to Reduce Ambient Light Noise in Visible-Light Communication using Solar Panels[END_REF] and the main expected benefits are to ideally obtain a self-supplied LiFi receiver [START_REF]Part 15.7: Short-Range Wireless Optical Communication Using Visible Light[END_REF], [START_REF] Haas | Wireless Data from Every Light Bulb[END_REF] operating outside and coupled with a large and flexible reception area. Remarkably, the use of organic solar cells as receiver for VLC has been demonstrated [START_REF] Zhang | Organic solar cells as high-speed data detectors for visible light communication[END_REF]. The latest performance of the commercially available ARMOR's OPV modules named ASCA® [START_REF] Armor Asca® | [END_REF], allows to further investigate the behavior of such technology. Therefore, the large scale availability of this third generation of PV technology is paving the way of future printed eco-friendly VLC systems. Some studies propose to use the PV cells left open-circuit (at least from a dynamic point of view) and one talks about open-circuit voltage (V oc ) operating mode. In this paper, we propose to provide performance results when the PV module N. Lorrière, M. Pasquinelli, L. Escoubas and J.-J. Simon are with Aix Marseille Univ, Univ Toulon, CNRS, IM2NP, Marseille, France (e-mail: nominoe.lorriere@im2np.fr; marcel.pasquinelli@univ-amu.fr; ludovic.escoubas@univ-amu.fr; jean-jacques.simon@univ-amu.fr).

G. Chabriel and J. Barrère are with Univ Toulon, Aix Marseille Univ, CNRS, IM2NP, La Garde, France (e-mail: gilles.chabriel@im2np.fr; jean.barrere@im2np.fr).

G. Pic and N. Vannieuwenhuyse are with ARMOR SAS company, La Chevrolière, France (e-mail: guillaume.pic@armor-group.com; nicolas.vannieuwenhuyse@armor-group.com).

is "shortcircuited" i.e when it is loaded by a transimpedance amplifier. We use the term: short-circuit current (I sc ) operating mode. Advantages of the I sc operating mode (vs. the V oc one) for higher data rate transmission are emphasized. Measurements of organic PV module efficiency under realistic outdoor conditions are provided. In particular, results about the impact of global and local irradiance [START_REF] Islim | The Impact of Solar Irradiance on Visible Light Communications[END_REF], [START_REF] Islim | An Investigation of the Solar Irradiance Effect on Visible Light Communications[END_REF] over the frequency response (gain and bandwidth) of a PV-based LiFi receiver are discussed for the two operating modes. Moreover, since a PV module offers a large detection surface, its behavior in response to a partial shadowing has to be explored. We show that shadowing influence is highly dependent of the operating mode and the shadow direction. II. MODELING I(V) characteristics of organic photovoltaic cells can be modeled in the same way as well-known inorganic cells. Indeed, the charge Density Of State (DOS) between HOMO and LUMO states is distributed exponentially [START_REF] Bässler | Charge Transport in Organic Semiconductors[END_REF], [START_REF] Clerc | Current status and challenges of the modeling of organic photodiodes and solar cells[END_REF] just as minority carrier on valence and conduction bands. Thus, in DC we consider that the relation between the output current I s and the output voltage V out is defined by the transcendental equation:

I s = I ph -I 0 exp V out + I s R s n kT q -1 - V out + I s R s R sh , (1)
with I ph the photogenerated current, I 0 the dark saturation current and kT q = V T the thermal voltage. The following figure Fig. 1 represents the classical linearized AC model corresponding to a forward biased solar cell. In this equivalent circuit, r d is the dynamic resistance, R sh the shunt resistance, C d the diffusion capacitance and R s is the lumped series resistance. Finally, i ph corresponds to the variations of the photocurrent.

A. V oc operating mode

When the output of the solar cell is left open (or loaded by a sufficiently high impedance), since i s ≈ 0, the effect of the series resistance R s can be neglected. Moreover, since a high DC current passes through the solar cell (I d ≈ I ph ), one can consider from

r d = nV T I ph , (2)
that the dynamic resistance is sufficiently low so that R sh can be neglected. Now, using

C d = τ r d , (3)
where τ represents the charge carrier lifetime, one obtains the following complex transfer function:

T oc (s) v out (s) i ph (s) v oc (s) i ph (s) = r d 1 + τ s . ( 4
)
It corresponds to a first-order low-pass filter where the static gain r d = nV T I ph is expected to decrease with the DC illumination I ph (and to increase with the temperature). At contrast, since the charge carrier lifetime τ generally decreases with the DC illumination, the -3dB cutoff frequency f c = 1 2πτ is expected to increase with the DC component of the photocurrent.

B. I sc operating mode

To use a solar cell in I sc operating mode, the cell has to feed the input of a transimpedance amplifier (TIA) as illustrated on the following figure Fig. 2, Since the OPA maintains a very low forward bias at the output of the solar cell, a very small portion of the DC current passes through the solar cell (say I d = εI ph , where ε 1). Consequently, the dynamic resistance r d becomes high and can then be neglected compared to R sh . Some mathematical derivations lead to the following three-order low-pass filter expression of the complex transfer function:

T sc (s) v out (s) i ph (s) (5) = A 0 Z i (s)Z F (s) (A 0 + 1 + τ 0 s)(R s + Z i (s)) + (1 + τ 0 s)Z F (s) where Z i (s) = 1 C d s R sh and Z F (s) = 1 C F s R F .
A rapid study shows that the static gain corresponds to the feedback resistance R F . It is then independent of the PV cell characteristics and therefore not linked to the DC injection level. Additionally, it can be shown that the maximum achievable cutoff frequency is 2) and ( 3), it comes

f c = 1 2πRsC d . From (
f c ∝ nV T R s εI ph τ . ( 6
)
In the above expression (6) of the maximum cutoff frequency, the DC current I ph appears in the denominator of the ratio and then tends to fall the cutoff frequency if it increases. Simultaneously, both the lumped series resistance R s and the charge carrier lifetime τ tend to decrease with the injection level I ph . Futher studies will be required to verify the variation of the cut-off frequency with DC illumination (the influence of temperature should also be considered).

III. EXPERIMENTAL SETUP

All the experimental frequency responses provided in this study were obtained using the OFDM (orthogonal frequencydivision multiplexing) testbench developed at IM2NP lab. and described in [START_REF] Lorrière | An OFDM Testbed for LiFi Performance Characterization of Photovoltaic Modules[END_REF]. The hardware consists, for the transmitter part, in a Philips Lumileds one High Power LED referenced as Luxeon Rebel Plus LX18-P140-Y. For the receiver part, we use an ASCA® organic PV module designed by ARMOR company (Fig. 3). It is composed of 9 cells all connected in series to exhibit an open circuit voltage of 7.28V and a short circuit current of 33.75mA at AM1.5 condition. The Transmitter/Receiver distance is set to obtain a 450 lux lighting at the LED maximum level illumination point. The data acquisition (AD conversion) is performed using a MDO3034 Tektronix oscilloscope and controlled by a MATLAB 2017b interface. For I sc mode measurements, a Stanford Research SR570 preamplifier (TIA) is added to convert the current provided by the module to a voltage suitable to the oscilloscope. The unipolar scheme chosen for the optical transmission is the classical DCO-OFDM (direct current biased optical OFDM) [START_REF] Dissanayake | Comparison of ACO-OFDM, DCO-OFDM and ADO-OFDM in IM/DD Systems[END_REF], [START_REF] Islim | Modulation Techniques for Li-Fi[END_REF]. Within each sent OFDM symbol, the K payload subcarriers k∆ f 1 ≤ k ≤ K, ∆ f = Fmax 2(K+1) , are phase-modulated by complex data extracted from a QPSK alphabet. Each symbol is extended by a cyclic prefix (CP) with a duration T cp . The cyclic prefix duration T cp has to be greater than the transient response duration of the channel in order to reach the steady-state regime. Here, T cp is chosen equal to the OFDM symbol duration i.e. 1 ∆ f . From K known transmitted data, averaging over M OFDM symbols and assuming a linear transmission channel, the orthogonal frequency modulation scheme allows to easily estimate, in one shot, the complex frequency response of the channel from the emitting LED (exhibiting a flat channel at considered frequencies) to the PV receiver for the whole subcarrier frequencies.

In this study, we get : K = 255, M = 975 and F max = 250kHz (resp. 50kHz) for the I sc mode (resp. for the V oc mode). For example, in the I sc mode, frequencies from

∆ f = Fmax 2(K+1) = 488.3Hz to K∆ f = 124.5kHz are scanned in 2 × M ∆ f ≈ 4s.
Obviously, the absolute value of the frequency response is not available. The provided gain frequency responses then have to be interpreted in a relative way from an experiment to another; the experimental conditions being maintained as closely as possible. Thus, since the two PV operating modes lead to two different experimental setups, the I sc and V oc measurements cannot be compared with each other.

IV. RESULTS

A. Influence of solar cell operating mode on the cut-off frequency 1) V oc operating mode: Fig. 4 provides the evolution of the gain frequency response for different values of the additional DC illumination level from 450 to 5500 lux.

As the DC level increases, it can be observed, as expected, an important diminution of the static gain T oc because r d ∝ 1 I ph (see Fig. 5 for a log-linear representation) and an increase of the frequency bandwidth f c ∝ 1 τ . As a consequence, in a V oc mode transmission context, the signal-to-noise ratio (SNR) is expected to degrade as the DC illumination level increases.

2) I sc operating mode: Fig. 6 shows the evolution of the gain frequency response for different values of the additional DC illumination level when LiFi communication is performed in I sc mode. Unlike Fig. 4, we can first observe that the static gain remains constant as expected by modeling. A slight increase of the cutoff frequency is also observed as the DC light level rises. The diminution of both the series resistance R s and the charge carrier lifetime τ could compensate the increase of the DC current I ph in the denominator of the expression of the cutoff frequency [START_REF]Part 15.7: Short-Range Wireless Optical Communication Using Visible Light[END_REF]. In addition to the non-sensitivity of the frequency response to DC illumination, it is important to note that without DC LED, the cutoff frequency on I sc mode (66.4 kHz) is 50 times higher than that on V oc mode (1.3 kHz). SNR measurements should help to confirm the net benefit of switching to I sc mode.

B. Shading effects

In this section, frequency responses are provided when the 110 mm × 52 mm area ASCA® OPV module is partially shaded in the lateral or longitudinal direction as illustrated on Fig. 7. It is well known [START_REF] Sullivan | Shadow effects on a series-parallel array of solar cells[END_REF], [START_REF] Hofer | Parametric analysis and systems design of dynamic photovoltaic shading modules[END_REF] that the two shading directions generally induce very different behaviors: when all the cells are partially shaded, it can be considered that it is just a PV module with a smaller active surface. When one or more cells are fully hidden, these non-illuminated cells act as impedance loads for the cells remained active. It is then clear that the shading location will differently impact the frequency responses according to the operating mode used:

• In short-circuit mode, since the amount of the received photocurrent i ph decreases, the output voltage R f i ph will decrease when all cells are more and more hidden (Fig. 8).

A high-pass filtering effect is observed (Fig. 9) when one or several cells are fully shadowed. This phenomenon is due to the shaded cells capacities in series acting as very high impedances for low frequency photocurrents. Moreover, the photocurrent i ph being generated by the only exposed cells, a general attenuation is observed. • In open-circuit mode, when all the cells are partially shadowed, the DC photocurrent I ph being weaker, the static gain (resistance r d ) tends to highly increase (non-linear behavior) and consequently tends to overcompensate the diminution of the dynamic photocurrent i ph . In the same time, the carrier lifetime increases (less available holes in the structure) and then diminishes the cutoff frequency. A higher voltage is then observed at low frequencies (Fig. 10).

When some cells are fully shadowed, most of the photocurrent is now absorbed by the non-illuminated cells. The modulus of these impedance loads (R sh C d ) decreasing with the frequency and decreasing when more and more cells are hidden, a diminution of both the static gain and the cutoff frequency is then observed (Fig. 11).

Generally, for shaded PV modules operating in the V oc mode, the output voltage decreases but not in the same degree as in the I sc mode. It can be explained by the increase of the static gain (r d ) witch is always constant (R f ) in the I sc mode.

V. CONCLUSION AND PERSPECTIVES

In this paper, we study in a first part the dynamic behavior of OPV modules as a function of the DC photocurrent for two operating modes: the V oc mode and the I sc mode. We show that in the V oc mode, the static gain diminishes with the DC illumination while the frequency bandwidth increases. A higher bandwidth is obtained in the I sc mode and both the static gain and the frequency bandwidth are a priori not impacted by an additional DC irradiation. This last result is encouraging for the use of PV modules as outdoor LiFi receiver and even more for mixed indoor/outdoor conditions. In a second part, we study the resilience of the two operating modes to shading effects. We have observed that the decrease of the static gain is more important in the I sc mode than in the V oc mode.

Finally, the I sc operating mode would allow to reach a better data rate than in the V oc mode but shading effects would impact the performance. Complementary SNR measurements are in progress to confirm such important results.

Fig. 1 .

 1 Fig. 1. PV cell AC low variations equivalent circuit

Fig. 2 .

 2 Fig. 2. Isc operating mode PV cell

Fig. 3 .

 3 Fig. 3. ASCA® OPV module

Fig. 4 .

 4 Fig. 4. Frequency response in Voc operating mode

Fig. 6 .Fig. 7 .

 67 Fig. 6. Frequency response in Isc operating mode

Fig. 8 .

 8 Fig. 8. Isc condition -Shading direction "% of all cells shaded"

Fig. 9 .Fig. 10 .

 910 Fig. 9. Isc condition -Shading direction "Number of cells shaded"

Fig. 11 .

 11 Fig. 11. Voc condition -Shading direction "Number of cells shaded"


Maxime Gobert
email: maxime.gobert@umons.ac.be

Franc ¸ois Vallée
email: francois.vallee@umons.ac.be

Jan Gmys
email: jan.gmys@umons.ac.be

Nouredine Melab
email: nouredine.melab@univ-lille.fr

Jean-Franc ¸ois Toubeau
email: jean-francois.toubeau@umons.ac.be

Daniel Tuyttens
email: daniel.tuyttens@umons.ac.be

Surrogate-Assisted Optimization for Multi-stage Optimal Scheduling of Virtual Power Plants

Keywords: Efficient Global Optimization, Two-stage optimization, power market, profit maximization, surrogate models

This paper presents a comparison between two surrogate-assisted optimization methods dealing with two-stage stochastic programming. The Efficient Global Optimization (EGO) framework is challenging a method coupling Genetic Algorithm (GA) and offline-learnt kriging model for the lower stage optimization. The objective is to prove the good behavior of bayesian optimization (and in particular EGO) applied to a real-world two-stage problem with strong dependencies between the stages. The problem consists in determining the optimal strategy of an electricity market player participating in reserve (first stage) as well as day-ahead energy and real-time markets (second stage). The decisions optimized at the first stage induce constraints on the second stage so that both stages can not be dissociated. One additional difficulty is the stochastic aspect due to uncertainties of several parameters (e.g. renewable energybased generation) that requires more computational power to be handled. Surrogate models are introduced to deal with that additional computational burden. Experiments show that the EGO-based approach gives better results than GA with offline kriging model using smaller budget.

INTRODUCTION

As industrial applications require more and more computational resources to achieve good levels of accuracy, we now face problems impossible to solve in a reasonable time. To handle this time constraint, many black-box-based optimization methods have been developed in the last decades [START_REF] Sacks | Design and analysis of computer experiments[END_REF]- [START_REF] Forrester | Recent advances in surrogate-based optimization[END_REF]. Black-box optimization is a global optimization process relying on surrogate models to substitute the so-called blackbox function which is a costly simulation that we do not need to know precisely. One of the state-of-the-art algorithms is Efficient Global Optimization (EGO), presented by Jones et al. in [START_REF] Jones | Efficient global optimization of expensive black-box functions[END_REF], on which is based this study. EGO is coupled with Gaussian Processes (GP), or kriging, as a surrogate model. Kriging seems to be a natural choice because of its particularity to provide an estimation of the prediction error, which is important to compute the Expected Improvement (EI). EGO with EI as infill criterion has proven to be a good choice for highly time-constrained optimization. It has been used in industrial applications such as the optimization of the shape of horn-loaded loudspeakers in [START_REF] Morgans | Ego shape optimization of horn-loaded loudspeakers[END_REF]. It also has been successfully applied to other research field like combinatorial problems in [START_REF] Zaefferer | Efficient global optimization for combinatorial problems[END_REF]. Literature points out the wide range of applications of bayesian optimization methods like EGO. Nevertheless, the application to multi-stage optimization is a very recent contribution. Multi-stage optimization has been rarely treated with surrogate models. In Sabater et al. [START_REF] Sabater | An efficient bi-level surrogate approach for optimizing shock control bumps under uncertainty[END_REF], a two-stage optimization is used for uncertainty quantification on the lower level. Uncertainty quantification is replaced by a surrogate model prediction and global optimization follows the maximization of EI scheme. In Islam et al. [START_REF] Islam | A surrogate assisted approach for single-objective bilevel optimization[END_REF], both stages are based on differential evolution (DE) and the inner stage is an optimization problem solved with surrogate model assistance for objective and constraint functions.

The proposed work builds on [START_REF] Toubeau | Medium-term multi-market optimization for virtual power plants: a stochastic based decision environment[END_REF] where the principal objective is to provide an efficient decision tool to maximize the total profit of an electricity market player participating in medium-term operating reserves as well as day-ahead and real-time energy markets. These platforms involve different embedded time horizons going from week-ahead (reserve market) to real-time. Consequently, the problem formulation takes the form of a two-stage algorithm where the first stage (reserve market decisions) is called Medium Term Optimization (MTO) and the second stage (decisions undertaken on the day-ahead and real-time markets) is called Short Term Optimization (STO). The STO is considered as a black-box taking as inputs the first stage decisions (reserve procurement) and as output a single real value: a variable profit, in euros. The distinctive characteristic here is that we do not solve the lower problem with a traditional optimization method. We only predict its optimal value using a metamodel trained with realizations of STO of which some first-stage variables (basically, the three existing electricity reserve products) are fixed at different (optimally-chosen) values at each iteration. Furthermore, the optimization tool must take into account the prediction uncertainty due to unexpected load or renewable production to contract energy reserves and be able to contribute to the grid safe state restoration in case of unexpected dizziness on the electricity transmission system. This is the reason why it is treated using a stochastic algorithm where several scenarios of possible short-term realizations are generated [START_REF] Toubeau | Deep learningbased multivariate probabilistic forecasting for short-term scheduling in power markets[END_REF].

This stochastic programming approach has two main challenges. The first one consists in modeling the uncertainty through a set of time-varying predictive scenarios that represent time trajectories of all uncertain variables. The second is to overcome the computational burden associated with the resulting formulation dealing with uncertainties. Both are treated in [START_REF] Toubeau | Medium-term multi-market optimization for virtual power plants: a stochastic based decision environment[END_REF]. This is why the use of a metamodel is mandatory to complete the optimization in the dedicated time (2 hours maximum due to electricity market rules).

To that end, in [START_REF] Toubeau | Medium-term multi-market optimization for virtual power plants: a stochastic based decision environment[END_REF], a kriging model is constructed as a preliminary task to be inserted in place of the STO. This type of substitution is often called offline learning and does not allow to improve the metamodel from experience (new simulations) and therefore, to control the evolution. In this paper, we investigate another type of surrogate-assisted optimization method to optimize the portfolio management. The approach is based on Bayesian optimization and more precisely on EGO [START_REF] Jones | Efficient global optimization of expensive black-box functions[END_REF] and qEGO variant from [START_REF] Ginsbourger | A Multi-points Criterion for Deterministic Parallel Global Optimization based on Gaussian Processes[END_REF]. We treat STO as a blackbox function and maximize the profit with a method able to improve the metamodel at each iteration. This kind of behavior is called online learning in opposition to offline learning. Experimentations are performed using a small but realistic case-study (Section III) in order to compare both approaches. The first one is constituted of an offline learnt kriging model from the SUMO Matlab toolbox [START_REF] Gorissen | A surrogate modeling and adaptive sampling toolbox for computer based design[END_REF] coupled to a Genetic Algorithm (GA) from Matlab Global Optimization Toolbox [START_REF]Matlab optimization toolbox[END_REF]. The second approach is the EGO-based one, relying on DiceOptim and DiceKriging R packages [START_REF] Roustant | Dicekriging, diceoptim: Two r packages for the analysis of computer experiments by krigingbased metamodeling and optimization[END_REF]. Experimental results show that the EGO-based approach finds improved solutions while requiring less calls to the simulator (namely: STO).

The paper is organized as follows. First, a short presentation of the two-stage stochastic optimization tool is provided in Section I. Theoretical aspects of the used methods are presented in Section II. Afterwards, in Section III we present the comparison of both approaches and display experimental results. Finally, some conclusions and perspectives for future works are drawn.

I. PROBLEM DESCRIPTION

A. Context and problem presentation 1) Context: Following the liberalization of the electricity sector, generation and transmission activities are fully decoupled. This context resulted in the development of aggregators, also referred to as Virtual Power Plants (VPPs) whose initial goal is to make some profit by helping the transmission system operator to maintain balance between generation and consumption (ensuring consequently a constant 50 Hz frequency within the grid). Practically, VPPs can be seen as single actors combining different generation, storage and load management technologies that are jointly operated within the objective to maximize their expected profit. In order to do so, VPPs can either participate in energy markets (selling/buying energy to other actors), or offer services to the transmission system operator to help him maintaining a safe and efficient system operation. Such services are typically contracted in mid-term (i.e. week or month ahead).

2) Multistage optimization: One can distinguish three levels in the optimization process: (1) the mid-term, (2) the day-ahead and (3) the real-time, each one is a market where decisions are taken. The combination of the day-ahead and real-time markets is also named short-term and is an optimization process itself. The decision tool is designed from the VPP perspective so that it maximizes its expected profit. The aggregator decision procedure involves two different (timedependent) stages. First, the optimization is done at the midterm level where decisions reserve procurement have to be taken. Practically, three products do exist depending on the required dynamics and are named R 1 (to be activated within 30 seconds to stabilize the frequency shift around 50 Hz), R 2 (to be fully activated within 7.5 and 15 minutes to help restoring the frequency to 50 Hz) and R 3 (longer-term reserve if frequency has not been set back to 50 Hz after use of R 1 and R 2 ). Let us call X u the decision vector of first stage variables. Second, based on that mid-term commitment, the contribution to day-ahead and real-time markets must be decided, it sets the second stage variables X l (e.g. activation of the reserves in real-time, unit commitment of generators). Since it is designed to couple two time horizons, it takes the form of a twostage optimization problem where optimization is made at each stage.

One could write the problem as follows:

max P t (X u , X l ) = max P f (X u ) + E [P v (X u , X l )] (1) s.t X l = argmax (E [P v (X u , X l )]) (2)
The X u variable is a real value triplet containing the volumes (R 1 , R 2 , R 3 ), and X l is a mixed-integer vector. P f and P v are the two contributions to the total profit P t , respectively corresponding to the fixed profit (linear combination of X u ) and the variable profit coming from the Mixed Integer Linear short-term Optimization problem (second stage) and related to X u .

B. Two-stage problem formulation

As pointed out above, the optimization of the two temporal horizons (mid-term and short-term) are intrinsically linked. Indeed, the resources allocated in mid-term decisions must remain available if they are to be activated in real-time. Fig. 1 illustrates both stages of the algorithm and their outputs P v and P f . Furthermore, the STO has to be carried out on a daily basis and some moderate risk attitude can thus be envisaged here. Two-stage scenario-based optimization is therefore favored over techniques such as robust or interval optimization techniques that are known to yield conservative (and thus suboptimal) solutions [START_REF] Bruninx | Improved modeling of unit commitment decisions under uncertainty[END_REF].

1) Mid-term Optimization: One of the major constraints is to maintain the balance on the power grid: at any time, the consumption must be equal to the production. Any failure to respect this balance can have disastrous effects such as congestion or blackout of the electricity system [START_REF] Toubeau | A new approach based on the experimental design method for the improvement of the operational efficiency in medium voltage distribution networks[END_REF]. The transmission system operator is responsible for this constraint and has to contract energy reserves to be triggered when needed. This situation instores an intricate relation between the electricity market and the reserve market so that it is highly preferable to tackle both problems simultaneously as a twostage optimization problem.

The first stage optimization (MTO) tries to determine the optimal participation to the services (reserves offered to the system operator). In practice, MTO decides the volume of reserves to be contracted. As already mentioned, the energy reserve is composed of three different types: primary, secondary and tertiary reserves(X u = (R 1 , R 2 , R 3 )). They are organized with respect to their temporal dynamics (response speed), with the primary product being the most reactive one. The total profit, the one we want to maximize, has two contributions: the first one called fixed profit depends linearly on the triplet fixed at mid-term P f (X u ), and the second is the result of the the non-linear second-stage optimization (STO) P v (X u , X l ).

2) Short Term optimization metamodeling: STO is an optimal scheduling problem and is itself a two-stage optimization problem (which makes the global problem a three-stage algorithm). It gathers the day-ahead and the real-time decisions. During the process, decisions concerning short-term market participation or reserve activations are taken regarding the midterm decisions. The decisions must consider constraints such as the physical limits of the various means of production, or the necessity to be balanced in real-time (i.e. what has been procured in previous market floors has to be actually delivered). The STO constitutes the non-linear part of the objective function: it gives the variable profit P v by optimizing the resource management in day-ahead and real-time markets. The short-term optimization is formulated as a Mixed Integer Linear Program (MILP) following the model from [START_REF] Pandžić | Virtual power plant mid-term dispatch optimization[END_REF]. It is subject to the three volumes of reserves set in the mid-term decision process.

II. SURROGATE-BASED OPTIMIZATION

Surrogate-based optimization (SBO) is used to tackle expensive-to-evaluate objective functions within an optimization process. The main idea is to momentarily sacrifice the model precision (using a surrogate model) in favour of execution speed. Using such a method allows to use global optimization methods that require a lot of evaluations -like population-based algorithms. Here, our attention is focused on a particular type of surrogate model, which is kriging.

A. Kriging

Kriging has been developed first in the field of geostatistics by the mining engineer Danie G. Krige in the early 1950s. The method has been formalized by Georges Matheron in the early 1960s [START_REF] Matheron | Principles of geostatistics[END_REF]. The kriging method consists in finding the best linear unbiased predictor (BLUP). It relies on a set of observations Y = (y (1) , . . . , y (n) ) at locations X = (x (1) , . . . , x (n) ) and their spatial correlation: how sensible is the target according to the location variation. The target Y is assumed to be the realization of a gaussian process N (µ(x), σ 2 (x)).

1) Kriging predictor: Let us consider the random variable Y (x) = µ(x) + Z(x) with µ(x) the trend, and Z(x) the residual function. The trend can be written as a linear combination of L basis functions and Z(x) as a centered gaussian process (GP) (as we removed the trend/mean). We end up with Y (x) = L l=0 a l φ l (x)+Z(x) and the BLUP at location x * is Ŷ (x * ) = n i=1 ω i Y (x (i) ) where Y (x (i) ) are the n observations. One (strong) hypothesis is that the centered gaussian process has to be stationary:

E[Z(x)] -E[Z(x + h)] = 0, ∀x, h ∈ R n and Var[Z(x) -Z(x + h)] = f (h) < ∞, ∀x, h ∈ R n .
This means that the mean is constant over the search space, and we have a spatial correlation between data. Intuitively, for two locations x 1 and x 2 , Z(x 1 ) must be close to Z(x 2 ) and only depends on h = u -v. This is the problematic part when dealing with optimization processes. We cannot be sure that for two executions of the optimizer with close parameters, we end up with close optimized solutions.

In the application targeted by this paper, the MTO only requires the variable profit resulting from STO. One of the challenges of this work is to substitute an optimization process (the STO) by a surrogate model. The problem lies in the nature of the optimization: for the given objective function of the STO P v : E ⊂ R n → F ⊂ R we are interested in finding an optimal point X l ∈ E which gives the optimum value of P v for a fixed X u . As the algorithm in question is a MILP one can not assume the continuity of ψ : X u → X l . But the fact is that we are searching for an optimal triplet

X * u = (R * 1 , R * 2 , R *
3 ) which maximizes the total profit P t . So that when considering the total profit computation P t (instead of P v ) as a black-box, it is reasonable to consider that for a neighborhood of X u , the profit remains in a neighborhood of P t (X u , X l ). Indeed, for close mid-term decisions, the optimal management of the remaining resources should earn similar variable profits and thus, similar total profits. Furthermore, the first approach consisting in building a kriging model as preliminary task shows good results [START_REF] Toubeau | Medium-term multi-market optimization for virtual power plants: a stochastic based decision environment[END_REF].

The BLUP must be unbiased and minimize the variance:

min Var[ Ŷ (x * ) -Y (x * )] such that E[ Ŷ (x * ) -Y (x * )] = 0.
The solution of the problem gives the universal kriging predictor which writes (from [START_REF] Roustant | Dicekriging, diceoptim: Two r packages for the analysis of computer experiments by krigingbased metamodeling and optimization[END_REF]):

m(x * ) = f (x * )ω T + c(x * ) T C -1 (Y -F ω) (3) σ(x * ) = C(x * , x * ) -c(x * ) T C -1 c(x * ) (4)
+ (f (x * ) T -c(x * ) T C -1 F ) T (F T C -1 F ) -1 (f (x * ) T -c(x * ) T C -1 F ) with ω = (F T C -1 F ) -1 F T C -1 Y
, and F = (f (x (1) ), . . . , f (x (n) )), f (x) being the vector of trend functions values at location x, C the covariance matrix and c(x * ) correlation vector. We will detail C and c(x * ) in the next section. Note that f (x * )ω T corresponds to the linear combination of basis functions: it is the trend value. The second part is related to the spatial correlation between data. More details on this formula can be found in [START_REF] Roustant | Dicekriging, diceoptim: Two r packages for the analysis of computer experiments by krigingbased metamodeling and optimization[END_REF].

2) Spatial correlation: Kriging is based on the spatial correlation between data which is represented by the covariance kernel. The covariance kernels used for this article are built upon the following expression from Roustant et al. [START_REF] Roustant | Dicekriging, diceoptim: Two r packages for the analysis of computer experiments by krigingbased metamodeling and optimization[END_REF]:

c(h) = C(u, v) = σ 2 d j=1 γ(h j , θ j ) ( 5
)
where

d is the dimension, h = u -v and σ = C i,i = C x (i) , x (i) = Var x (i) .
Hence, from the previous statement, we have for equations (3): C = C x (i) , x (j) i,j , i, j = 1, . . . , n and c(x * ) = C x * , x (i) i , i = 1, . . . , n. And θ is a vector of hyperparameters that must be fitted. This is often done by numerical Maximum Likelihood Estimation (MLE), which is an optimization process.

We usually choose a covariance kernel among well known forms. Indeed, it must satisfy some constraints to get the resulting covariance matrix positive definite. One can find those conditions in N. Cressie's book [START_REF] Cressie | Statistics for Spatial Data[END_REF]. Verifying such an assumption is not an easy task and this is why we usually choose a kernel among a few forms. Among often used kernel functions, we can find: exponential, gaussian or matérn kernels. Their analytical form can be found in [START_REF] Roustant | Dicekriging, diceoptim: Two r packages for the analysis of computer experiments by krigingbased metamodeling and optimization[END_REF].

This model is often presented as the DACE model -Design and Analysis of Computer Experiments, popularised by Sacks et al. in [START_REF] Sacks | Design and analysis of computer experiments[END_REF].

B. Efficient Global Optimization (EGO)

1) EGO algorithm: EGO is an online learning algorithm, meaning that the metamodel is improved with time. But even more, it actively searches for interesting areas to explore and selects the most promising point. This algorithm is based on a metric, called acquisition function, which selects the point that most likely improves the predictor. We can mention the Probability of Improvement (PI) and the Expected Improvement (EI) [START_REF] Schonlau | Computer experiments and global optimization[END_REF]. They both rely on the uncertainty of the predictor, which is available in the kriging model.

The EGO algorithm works through several stages repeated sequentially.

• First, let us consider a sample and fit the model (kriging model in our study) to the data. We need a sample that represents the search space fairly well. Jones [START_REF] Jones | Efficient global optimization of expensive black-box functions[END_REF] advises to start with ≈ 10d points from Latin Hypercube Sampling (LHS), where d is the problem dimension. Another recommandation comes from Sóbester et al. [START_REF] Sóbester | On the design of optimization strategies based on global response surface approximation models[END_REF] who suggest to initiate the process with 35% of the budget.

• The second step is to maximize the chosen acquisition function: EI, defined by Equation [START_REF] Morgans | Ego shape optimization of horn-loaded loudspeakers[END_REF].

EI(x) =(y opt -y pred (x))Φ y opt -y pred (x) σ(x) (6) +σ(x)φ y opt -y pred (x) σ(x)
where: Φ is the probability density function (pdf) of N (0, 1) and φ is the corresponding cumulative density function (cdf). EI is high when the prediction is close to the current known optimum or when the uncertainty is large. Conversely, if the prediction is reliable or not attractive, EI will be low.

There exists a few methods for that purpose, we can mention the Branch-and-Bound algorithm from Jones et al. [START_REF] Jones | Efficient global optimization of expensive black-box functions[END_REF] and the genoud and BFGS algorithms from the R package DiceOptim [START_REF] Roustant | Dicekriging, diceoptim: Two r packages for the analysis of computer experiments by krigingbased metamodeling and optimization[END_REF]. If the maximization of EI is able to provide a new point (not too close to an already known area), this latter is evaluated and the model updated.

• Once the point maximizing EI is found, it is evaluated with the real simulator (STO). This new known point is added to the sample and the model is updated.

The process is repeated until the stopping criterion is met. Stopping criterion can be for example a given budget, a maximum time reached, or a condition on the improvement. Note that the update of the model as well as the maximization of EI can become costly when the sample becomes important. Algorithm 1 summarizes the process through the pseudo-code and Fig. 2 gives an illustration of one EGO cycle.

2) qEGO adaptation: Another version of the EGO algorithm has been developed by Ginsbourger et al. in [START_REF] Ginsbourger | A Multi-points Criterion for Deterministic Parallel Global Optimization based on Gaussian Processes[END_REF]. The structure of the algorithm remains the same but the selection method of new points is modified. Several points may be selected (q points) to be added to the sample. We present here only the practical method to get a multi-point expected improvement. The interested reader can find the theoretical description of this method in Ginsbourger et al. [START_REF] Ginsbourger | A Multi-points Criterion for Deterministic Parallel Global Optimization based on Gaussian Processes[END_REF].

Algorithm 1 EGO pseudo-code Create initial design X = [x (1) , . . . , x (n) ], Evaluate Y = f (X), Fit the kriging model, while stoping criterion = false do x (n+1) = argmax D (EI(x)) Update the kriging model adding (x (n+1) , y (n+1) ) end while One convenient way to compute q-points EI is to run sequentially q EGO steps without the costly evaluation. This leads to the following heuristics to attribute a value to the selected point.

a) qEI heuristics:

• Kriging Believer (KB) method consists in trusting the metamodel prediction. We run EGO, select one point and attribute the predicted value as it was the answer from the simulator. We learn a new model based on the n known points to which we add the best candidate, we end up with n + 1 points. This sequence is reapeated until we get the q points and only then the q elements are evaluated.

Believing in the kriging accuracy can misguide the optimization process because of the current quality of the model. Once again, the theoretical uncertainty on the kriging predictor can be used to counter that drawback, which leads to the Constant Liar heuristic.

• Constant Liar (CL) method is very similar to the previous one, the only difference is that we use the confidence interval given by the kriging surrogate model. Instead of using the prediction, we can decide if we use the upper or the lower prediction bound. Let us denote those two methods CL min and CL max .

Say we are searching for a minimum, CL min has an attractive effect on the exploration process and then we favor intensification. Conversely CL max has the opposite effect and then we favor exploration. One could imagine randomly use one or another version of CL heuristic to balance exploration and exploitation.

III. QEGO APPLIED TO OPTIMAL SCHEDULING OF VPP

A fair comparison requires to take into account both execution time and simulation budget. Indeed, it is not appropriate to account only for evaluation time, considering that the simulation time is strongly impacted by the aggregator portfolio and we are dealing with a small actor (around 30 seconds per evaluation) so that 2 hours is more than enough to deal with such a case. We also observed that most of the computation time is spent for the simulations and that the extra cost of learning a new kriging model at each cycle is not prohibitive to consider the simulation budget as a fair comparison metric. This is particularly true for longer simulations.

A. Settings of the testing environment 1) Hardware and software: Experimentations are run on a single core of an Intel(R) Core(TM) i7-7500U CPU clocked at 2.70GHz. We choose the R package DiceOptim [START_REF] Roustant | Dicekriging, diceoptim: Two r packages for the analysis of computer experiments by krigingbased metamodeling and optimization[END_REF] because it is user friendly and sufficiently complete. Indeed, DiceOptim is based on DiceKriging which contains several kernel functions, and various kriging forms for preliminary studies. Furthermore, in addition to classical EGO algorithm, DiceOptim contains functions related to qEGO. The use of qEGO, even in a sequential mode can be advantageous because of the CL heuristic presented above. It can be used as an additional parameter to balance exploitation and exploration. The simulator is written in Matlab so we interface R and Matlab using the RMatlab package from R [START_REF] Henrik | Read and write mat files and call matlab from within r[END_REF]. The second compared method is implemented with Matlab and the Global Optimization toolbox [START_REF]Matlab optimization toolbox[END_REF] for the GA. The surrogate is constructed thanks to the SUMO toolbox [START_REF] Gorissen | A surrogate modeling and adaptive sampling toolbox for computer based design[END_REF].

2) Description of the VPP test case: The VPP is composed of 3 Conventional Power Plants (CPP), with a maximum output power of respectively 130, 80 and 55 MW as well as 2 Pump Storage Units (PSH), both characterized by an output power of 15MW (but with energy limitations of respectively 25 and 75 MWh). Then, there are also renewable generation, i.e. wind farms and domestic rooftop photovoltaic (PV) installations totalizing 120MW of power.

The mid-term uncertainty is addressed by defining statistically representative days of wind and solar generation as well as total consumption within the portfolio. Afterwards, the day-ahead scenarios are generated. The VPP is paid for both the availability of reserves (capacity) and their actual activation (energy). The prices for the reserve capacity are fixed to 16e/MW for R 1 , 4e/MW for R 2 , and 1e/MW for R 3 . The activation prices reflect the technology-specific operation costs. The portfolio is arbitrarily created based on real data from the Belgian system in order to represent a typical actor. The generation and consumption patterns are realized regarding aggregated data for a typical month of July.

B. Technical aspects and initial settings of the methods

Both algorithms have parameters (or hyperparameters) that must be tuned for each algorithm.

1) Genetic algorithm with Gaussian Processes: The inner part (STO) is replaced by a kriging model trained with a classical Latin Hypercube Sampling to which the 8 corner points are added for a total of 48 points. The mid-term optimization is the part of the algorithm which determines the optimal triplet (R * 1 , R * 2 , R * 3 ) that realizes the maximum total profit. MTO is performed by a GA from the Matlab Global Optimization toolbox [START_REF]Matlab optimization toolbox[END_REF]. Each individual of the GA is encoded by (R 1 , R 2 , R 3 ). The GA initializes the population with a uniform sampling. The selection operator is stochastic uniform. It attributes a section of a line of length proportional to the individual score and randomly picks a number to select the interval. The crossover function is a scatter function that generates a random binary operator and selects attributes from parent 1 or 2 according to the corresponding binary element, its probability is 0.8. Mutation is a small perturbation of each entry of the individual taken from a centered gaussian distribution. The gaussian standard deviation is set to

σ k = σ k-1 (1 -k N ),
where k is the current generation, N the maximum number of generations, and σ 0 = 1. Its probability of appearance is 0.01. The number of generations is set to 100d, which makes 300 generations maximum. Each generation is constituted of 50 individuals. The previous parameters constitute the default choices of the toolbox.

2) Efficient Global Optimization: According to preliminary experiments we decided to start with a relatively small sample containing 15 points. This roughly corresponds to 30% of the simulation budget (which is 47). 15 points are enough to represent fairly the landscape and to realize the first kriging model. The initial sampling is done by a LHS from the R package DiceDesign [START_REF] Dupuy | DiceDesign and DiceEval: Two R packages for design and analysis of computer experiments[END_REF].

The kriging model uses a linear trend and a Matérn covariance kernel. Linear trend is chosen because of the linear part of the total profit, and the kernel choice is based on literature. Likelihood maximization is done with the genoud algorithm because it gives better results than the other algorithm proposed in the package, i.e. the BFGS algorithm. For larger problems which necessitate more sample points, BFGS might be favored because it is faster. Concerning the optimization process, we select the CL heuristic for the qEI points selection. This method presents interesting behavior for the exploration balance [START_REF] Ginsbourger | A Multi-points Criterion for Deterministic Parallel Global Optimization based on Gaussian Processes[END_REF]. The maximization of EI is done by the genoud algorithm which is the defaut parameter in DiceOptim [START_REF] Roustant | Dicekriging, diceoptim: Two r packages for the analysis of computer experiments by krigingbased metamodeling and optimization[END_REF]. A brief analysis on benchmark functions revealed that even in sequential, qEGO can favor diversification in the search and avoid to remain stuck in local maxima. To confirm that, we investigate qEGO for the treated problem.

C. Experimental results

1) Investigating qEGO: The goal of this paragraph is to evaluate the impact of the parameter q on the solution quality. Fig. 3 presents the behavior of EGO and qEGO with q = 2, 4, or 8. The plots show the best known target (i.e. the total expected profit) relatively to the size of the sample, this means that for q = 2 we have 16 cycles, q = 4 allows 8 cycles and q = 8 only 4, against 32 for traditional EGO. Regarding this, it is not surprising that for equivalent sample size, the best known target is not as good for q = 8 which has 4 times less surrogate updates than for q = 2. But it is also a good point to get a wider exploration which is indicated by a higher standard deviation (for equivalent sample size).

The analysis presented in Fig. 3 has been realized with the same initial sample for every search. It shows that q = 2 seems to be a good compromise to explore wider without delaying the convergence. Indeed, a larger box in the boxplot reveals a larger variance and so a wider exploration. If we convert the sample size parameter in number of cycles we observe that for less cycles we also achieve good results. In particular, we reach a good expected profit with q = 8 with only 4 cycles. Let us remind that a cycle covers a single surrogate model update, and q simulations.

2) Comparison of the two approaches: qEGO against GA + kriging: All the executions of the SUMO approach took around 25 minutes each, and for the same simulation budget the qEGO algorithm took around 30 minutes. Even if the number of calls to the simulator is the same (47 in total, including initial training set), we observe a few minutes difference due to the update of the kriging model at each qEGO cycle and to the maximization of EI.

Fig. 4 presents a boxplot realized with 10 executions of the qEGO algorithm running with different initial samples. We first observe a fast improvement of the target (P t ) even though the standard deviation remains important. We also can see that after a few cycles (q = 2 points added per cycle), the variance shrinks and does not vary much after 10 cycles, which can be interpreted as a convergence to an optimum. The algorithm does not improve significantly the optimal target with more cycles. The qEGO-based approach reveals to be slightly more time consuming for equal budgets. Nevertheless, it is important to notice that the optimal solution area is found by qEGO within the 10 first cycles (cf. Fig. 4) and the algorithm would have given a satisfactory result (better than GA-based approach) even with truncated execution, this result is evidenced by Fig. 5.

Table I presents the expected optimal profit realized by the aggregator. The first column shows qEGO results, which are truly evaluated by the simulator. The second column displays GA results within two sub-columns. The GA-based method runs with the surrogate model that returns predicted solutions. Optimal predicted solutions are truly evaluated a posteriori and presented in the real column. One can notice a substantial difference between the real and predicted values.

If we compare the final target of the qEGO-based approach with the one from GA-based approach on Table I, we can see that on average the daily expected total profit obtained by the qEGO approach is 244e greater than the one obtained by the GA-based approach. This corresponds to an average annual gain of 89000e. We must also highlight the very small standard deviation of the qEGO results that show the robustness of the algorithm. Indeed the GA-based approach gives very different results from an execution to another and often overestimates the optimal value that results in a suboptimal choice. When considering the extreme case realizations, we find an annual differences ranging from 50 000e to 155 000e. Let us remind that this is a relatively small actor so that one can expect much bigger benefits with major actors. Fig. 5 presents the results of both methods in the form of a boxplot. One can see that for any qEGO execution, we obtain a better final expected profit (qEGO18 box) than for any execution of the GA-based approach (GA SUMO box). Furthermore, the small area of the box reveals a small standard deviation and so a stable result. Other presented boxes correspond to prematurely stopped executions of qEGO: qEGOa represents the expected total profit after a cycles. It confirms that the qEGO algorithm is more efficient in finding an optimal expected profit, and qEGO would have given good results (still better than GA+SUMO) with an early stop.

CONCLUSION AND FUTURE WORK

To overcome issues in the electrical transmission network, the system operator must contract reserves to one or more aggregators (also referred to as VPPs). Once it is settled, the VPP has to uphold the contract and reserves that are contracted in mid-term must remain available in case of need in real-time. Non-respect of the contract implies financial penalties for the VPP. Thus, omission of the dependencies between different time horizons can lead to non-feasible or non-tractable solutions. To deal with this issue, a two-stage stochastic optimization algorithm has been developed. The computational cost of the lower level requires the use of surrogate models to be able to treat the optimization within the time limit of 2 hours. The first approach presented here and previously developed in [START_REF] Toubeau | Medium-term multi-market optimization for virtual power plants: a stochastic based decision environment[END_REF] consists in coupling a GA to an offline learnt kriging model which raplaces the STO. The second approach is based on the EGO framework. This is an online learning method that presents the advantage of having a population constituted only of simulated individuals.

The metamodel starts with a restricted budget and improves with new simulations that are given by the maximization of the EI metric. Another advantage is that this algorithm is able to provide a solution even if the execution is interrupted (e.g. due to time constraints).

Considering a small, but realistic case study, the presented experimental results show that EGO and qEGO outperform the GA-based approach in terms of solution quality and execution time. Indeed, the VPP's expected profit is always larger with the qEGO-based approach. Furthermore, the optimization result is very stable between executions (very low variance). The expected extra profit for the VPP is evaluated up to 155 000e (+3.4% compared to GA+kriging) for a year when using the qEGO-based approach for the studied testcase. Good results of the proposed optimizer coupled with surrogate models makes realistic the management of bigger problems. As mentioned, the aggregator for this study might be a real case, it remains a relatively small actor. Considering bigger actors (with bigger portfolios) and the multiplication of renewable energy sources will increase the computational burden. Therefore, combining multi-level massively parallel computing with surrogate-assisted optimization (qEGO) will be a viable and important perspective of this work.

Fig. 1 .

 1 Fig. 1. Summary of the two-stage problem description

Fig. 2 .

 2 Fig. 2. One step of EGO illustrated

Fig. 3 .Fig. 4 .

 34 Fig.3. Influence of the q parameter for qEGO algorithm

Fig. 5 .

 5 Fig. 5. Boxplot of 10 realizations for the two approaches, with truncated results of qEGO (q=2)


Marion Princaud

Cyril Aymonier

Anne Loppinet-Serani

Nicolas Perry

Guido Sonnemann

Environmental Feasibility of the Recycling of Carbon Fibers from CFRPs by Solvolysis Using Supercritical Water

Keywords: Life cycle assessment (LCA), Supercritical water, Solvolysis, Recycling, Composites, CFRP

Originally developed for high-tech applications in the aeronautic and aerospace industry, carbon/epoxy composites have been increasingly used in the automotive, leisure, and sports industries for several years. Nevertheless, the carbon reinforcement is an expensive constituent, and it has been recently shown that it is also the most environmentally impacting in a composite part manufacturing. Recycling these materials (even restricted to the reinforcement recovery) could lead to economic and environmental benefits, while satisfying legislative end-oflife requirements. The solvolysis of the matrix by water under supercritical conditions is an efficient solution to recover the carbon fiber reinforcement with mechanical properties closed to the ones of virgin fibers. This paper aims at demonstrating the environmental feasibility of the recycling of carbon fiber/thermoset matrix composites by solvolysis of the matrix in supercritical water. This demonstration is based on life cycle assessment that evaluates benefits and environmental challenges of this recycling loop.

■ INTRODUCTION

Carbon fiber-reinforced plastics (CFRPs), or thermoset matrix composites, were originally developed for high-tech applications in the aeronautic and aerospace industry. For several years now, these materials have also been increasingly used in the automotive, leisure, and sports industries. In many applications in these sectors, one may seek aesthetic criterions or a simple feeling of high technology, more than highly technical properties. Thus, constituents' characteristics, and specifically reinforcements, are considered as a secondary matter and may be overemphasized regarding the function of the product. This is particularly true for nonstructural decorative parts (e.g., with a carbon look finish), for which the reinforcement is the most expensive constituent, and where glass fibers, much more less expensive, cannot be used. [START_REF] Pimenta | Recycling carbon fibre reinforced polymers for structural applications: Technology review and market outlook[END_REF][START_REF] Mcdonough | Applying the principles of green engineering to cradle-to-cradle design[END_REF] Today, there is no, or a limited, deposit (or very few) of carbon fibers from airplanes at the end of life because airplanes integrating such materials are only currently being built and will become waste later. In the future, the expected amount will grow year after year. Therefore, the question is this: Could carbon fibers recycled from airplanes (or from production waste from aircraft and automotive production) substitute mechanically for the majority of carbon fibers currently used in the automotive, leisure, and sports industries, considering that the recycling can be done in a cost-effective way and that the aeronautic industry will not use recycled fibers? Subsequent questions are these: How can carbon fiber-reinforced plastics be recycled? Is the recycling environmentally more sustainable than the production of virgin carbon fibers?

One of the first uses of the supercritical fluid technology in the field of recycling was applied to polymers. This technique has been developed extensively in Japan since 1995 and has been reviewed many times. [START_REF] Goto | Chemical recycling of plastics using sub-and supercritical fluids[END_REF][START_REF] Adschiri | Green materials synthesis with supercritical water[END_REF][START_REF] Aymonier | Supercritical water for environmental technologies[END_REF] Beyond plastics recycling, solvolysis in near-and supercritical fluids of thermosetting resins (phenol and epoxy resins) has attracted a great interest among the scientific community to recover materials like carbon fibers with a high added value in the past few years. To date, few studies have been carried out on the chemical recycling of these waste composites with near-and supercritical solvolysis technology. [START_REF] Morin; Loppinet-Serani | Near suband supercritical solvolysis of carbon fibres reinforced polymers (CFRPs) for recycling carbon fibers as a valuable resource: State of the art[END_REF][START_REF] Sugeta | Decomposition of fiber reinforced plastics using fluid at high temperature and pressure[END_REF][START_REF] Kamimura; Yamada | Effective depolymerization waste FRPs by treatment with DMAP and supercritical alcohol[END_REF][START_REF] Okajima; Yamada | Decomposition of epoxy resin and recycling of CFRP with sub-and supercritical water[END_REF][START_REF] Pinero-Hernanz | Chemical recycling of carbon fibre reinforced composites in nearcritical and supercritical water[END_REF][START_REF] Jiang | Characterisation of carbon fibres recycled from carbon fibre/epoxy resin composites using supercritical n-propanol[END_REF][START_REF] Pinero-Hernanz | Chemical recycling of carbon fibre composites using alcohols under subcritical and supercritical conditions[END_REF][START_REF] Liu | Recycling of carbon fibre reinforced composites using water in subcritical conditions[END_REF][START_REF] Saito | Science and Technology of Supercritical Fluids[END_REF] Compared to other recycling processes (mechanical recycling processes, pyrolysis, fluidized bed processes, low temperature solvolysis processes), near-and supercritical solvolysis has the huge advantage that clean carbon fibers are recovered with similar mechanical properties to pristine fibers. [START_REF] Morin; Loppinet-Serani | Near suband supercritical solvolysis of carbon fibres reinforced polymers (CFRPs) for recycling carbon fibers as a valuable resource: State of the art[END_REF] Moreover, these undamaged fibers are obtained at relatively low temperature, without using organic solvents or concentrated acids.

Near-and supercritical water and alcohols were mainly processed as solvolysis media. In fact, near-and supercritical water or alcohols play the role of solvent and reagent for the depolymerization of condensation polymers by solvolysis into their monomers in a fast and selective way; it will be a hydrolysis reaction with water and an alcoholysis reaction with alcohols. Condensation polymers are constituted with ether, ester, or acid amide linkages, which can be broken by hydrolysis or alcoholysis. The example of polyethylene terephthalate (PET) bottle recycling is significant in term of quantity but also of development of supercritical fluid-based recycling technologies. PET can be hydrolyzed in terephthalic acid (TPA), its monomers, and ethylene glycol in sub-and supercritical water. [START_REF]Recycling Thermoset Composites of the SST. 7th PCRD[END_REF] Composite plastics such as glass and carbon fiberreinforced plastics can be decomposed into monomers and fiber materials. Some years ago, the successful hydrolysis of an isolated epoxy resin in sub-and supercritical water has been already carried out. [START_REF] Fromonteil | Hydrolysis and oxidation of an epoxy resin in sub-and supercritical water[END_REF] The solvolysis of composite materials using near-and supercritical fluids, especially water and alcohol, was recently reviewed by our group. [START_REF] Morin; Loppinet-Serani | Near suband supercritical solvolysis of carbon fibres reinforced polymers (CFRPs) for recycling carbon fibers as a valuable resource: State of the art[END_REF][START_REF] Loppinet-Serani | Hydrolysis in Near-and Supercritical Water for Biomass Conversion and Material Recycling[END_REF] This way is very efficient in a technological point of view, but what is about sustainability?

In this paper, the results from an initial life cycle assessment of the supercritical fluid technology applied to CFRPs recycling is proposed for the first time in order to position recycled carbon fibers in its market between virgin carbon fibers and glass fibers with regard to sustainability and cost considerations.

■ MATERIALS AND METHODS

Environmental Assessment. Data for CFRPs Composites Manufacture. Duflou et al. [START_REF] Duflou | Environmental impact analysis of composite use in car manufacturing[END_REF] have some life cycle assessment (LCA)-based information on the environmental impacts due to petrochemical manufacturing of composite parts for vehicles as an alternative to steel, for lightening the vehicle, and for reducing life cycle air emissions beyond the benefits of plug-in vehicles. [START_REF] Michalek | Valuation of plug-in vehicle life-cycle air emissions and oil displacement benefits[END_REF] In a conventional car, the use phase has the greatest environmental impact due to high fuel consumption (directly related to the mass of the vehicle). In its lighter alternative version, it is the manufacturing phase that could become predominant. [START_REF] Michalek | Valuation of plug-in vehicle life-cycle air emissions and oil displacement benefits[END_REF] This is due to the carbon fiber manufacturing (see our analysis in Figure 1) based on data from Duflou et al. [START_REF] Duflou | Environmental impact analysis of composite use in car manufacturing[END_REF] and recalculated relative to the mass of the chosen product, i.e., 1 kg of carbon fiber. Furthermore, the main source of impact for these carbon fibers is due to the use of fossil fuel that has an important carbon foot print. [START_REF] Grice | Life cycle carbon footprint of re-refined versus base oil that is not re-refined[END_REF] Hence, it might be of real interest from a sustainability point view to propose recycled fibers as a way forward to limit the environmental impacts of the composite parts of light cars.

Due to the fact that the carbon reinforcement is the most impacting constituent in a carbon/epoxy composite's elaboration process (Figure 1), [START_REF] Duflou | Environmental impact analysis of composite use in car manufacturing[END_REF] recycling end-of-life composites (even restricted to the reinforcement recovery) could lead to reduce some anthropogenic impacts by decreasing the use of first-generation raw materials (mainly petroleum) for their production. Besides, it would help design engineers to balance energy efficiency and cost, by opening new opportunities for developing second-generation composites first dedicated to the manufacture of medium or low loaded parts. Lastly, recycled carbon fabric could widen the range of reinforcements on the marketplace between first-generation carbon and glass fibers.

All this has to be done in line with European directives that already force industries to improve their products' recyclability (e.g., in automotive industry 22 ). However, making feasible this new recycling sector requires overcoming users' reluctances by ensuring the secondgeneration semi-product's validity from economic and environmental aspects. Therefore, we carried out a life cycle assessment (LCA) in which the resource efficiency and potential environmental challenges of the carbon/epoxy composites' recycling process are analyzed.

Life Cycle Assessment: Goal and Scope. Every stage of the life cycle of the composite part has to be modeled in the LCA, from its manufacture to its end-of-life treatment, following the usual steps defined by the ISO 14040 standards. 23 These ISO standards define LCA as the following: "Compilation and evaluation of the inputs, outputs and the potential environmental impacts of a product system throughout its life cycle". LCA is the only method that assesses the environmental impacts of a product or activity over its entire life cycle. It is a holistic approach that takes into account the extraction and treatment of raw materials, product manufacturing, transport and distribution, and product use and end-of-life. LCA is structured in the following phases: (a) goal and scope definition, (b) life cycle inventory, (c) life cycle impact assessmen, and (d) interpretation. Life cycle impact assessment assigns life cycle inventory results to impact categories like climate change and ionizing radiation; the environmental profile consisting of the indicator results for the impact categories selected provides information on the environmental issues associated with the inputs and outputs of the product system under study.

As previously mentioned, we focus on carbon/epoxy composites. The resin is an epoxy one. The carbon fibers were furnished by industry partners; therefore, we do not have any information about their precise nature. The deposit of materials to be recycled consists possibly in end-of-life aeronautic parts but, most likely to date, in composite offcuts. The composite part chosen for the LCA is assumed to be processed in Europe with Japanese carbon reinforcement. Its mass is supposed to be 1 kg. Thus, we aim at studying the interest of recycling such materials more generally, such as the environmental feasibility of the recycling process.

Life Cycle Inventory. The following analysis is based on Duflou's data, [START_REF] Duflou | Environmental impact analysis of composite use in car manufacturing[END_REF] which assessed the manufacturing of composite semi-structural panels in automotive industry. All of these data have been recalculated relative to the mass of the chosen product (i.e., 1 kg). In our case study, the use phase is not taken into account. Indeed, to the best of our knowledge, the only input data that can be taken into account concern transport operations. Like so, as rather classically, the present simulation shows that this factor did not contribute much to the overall impacts (less than 5%).

Regarding the product's end-of-life, two scenarios have been modeled: The first one consists of burying the composite part, which is what is currently done, and represents the reality for actual composites at their end-of-life. The second one consists of the recovery of the carbon reinforcement. We focus on the recycling process by solvolysis described in Figure 2). We consider (i) an aqueous solvolysis of the matrix by water under supercritical conditions (temperature around 400 °C and pressure about 25 MPa) and (ii) a hydrothermal oxidation of the effluent to clear matrix components from water at the end of the solvolysis process.

This technology allows the fiber to be recovered. Therefore, it is a real (but partial) recycling and not a simple material valorization. [START_REF] Morin; Loppinet-Serani | Near suband supercritical solvolysis of carbon fibres reinforced polymers (CFRPs) for recycling carbon fibers as a valuable resource: State of the art[END_REF] Lastly, the process uses energy, water, and oxygen, and only emits water and carbon dioxide.

Lastly, the research team from the Mechanics Institute of Bordeaux has developed a prototype for packaging these second-generation fibers in an attractive form for users (i.e., designers). Data matching the remanufacturing stage have not been taken into account yet in this very first LCA. However, this energy input is assumed to be very weak compared to those involved in the first-generation reinforcement process. As a consequence, the life cycle only loops after the manufacturing of the first-generation carbon reinforcement, with no specific additional remanufacturing.

Life Cycle Assessment: Software, Database, and Method. The LCA is carried out with the SimaPro software (v.7), [START_REF] Goedkoop | SimaPro 7 -Tutorial[END_REF] Eco Invent database (v.2), [START_REF] Frischknecht | The ecoinvent database system: A comprehensive web-based LCA database[END_REF] and ReCiPe Midpoint (H) method. [START_REF] Alcimed Report | Evaluation du gisement de fibres de carbone àrecycler et de leur valorisation pour le consortium RECCO[END_REF] As previously mentioned, in the recycling stage, the avoided material is the reinforcement. In other words, the production of a new raw material with nonrenewable resources (i.e., first-generation carbon reinforcement) is avoided. MPa. Water in the reactor was found to inhibit the coking and enhance the decomposition of the resin compared with the case of pyrolysis. As a result, clean carbon fiber was recovered, and the resin was decomposed and removed from the carbon fiber. [START_REF] Okajima; Yamada | Decomposition of epoxy resin and recycling of CFRP with sub-and supercritical water[END_REF] It can be pointed out that this solvolysis process is able to treat all types of composites, no matter their surface quality, geometry, size, density, etc. The only constraint is the reactor geometry.

■ RESULTS

Recycled Carbon Fibers

Morin et al. have also performed the recycling of carbon fibers from carbon fiber-reinforced composites in a semicontinuous flow reactor. Experiments were carried out at a temperature around the critical temperature of water for a reaction time of about 30 min. The process has been optimized in order to improve the solvolysis rate of the resin without the degradation of the mechanical properties of the fibers. Water or alcohols can be used as the solvolysis medium. They are different in terms of energy consumption because the critical coordinates of alcohols are generally lower than those of water. Therefore, recycling of CFRPs using an alcoholysis process could require less energy, but the hydrolysis process is safer and greener. In this study, water was used as solvent for the recycling of carbon fibers from CFRPs. The epoxy resin was completely decomposed into lower molecular weight organic compounds. Recovered carbon fibers were characterized using thermogravimetric analysis (TGA) to determine the amount of resin removed by the process, scanning electron microscopy to observe the fibers, and single fiber tensile tests to evaluate the mechanical properties of the recycled fibers. Recycled carbon fibers from CFRPs are clean (Figure 3). All the resin was removed according to the TGA results. Furthermore, the recycled carbon fibers present good mechanical properties; a tensile loss close to the one of virgin fibers is obtained. [START_REF] Morin; Loppinet-Serani | Near suband supercritical solvolysis of carbon fibres reinforced polymers (CFRPs) for recycling carbon fibers as a valuable resource: State of the art[END_REF] The final liquid phase is also analyzed by gas chromatography. The monomers of the initial resins have been identified.

Environmental Evaluation. The LCA of a 1 kg composite part that takes into account the recycling of the reinforcement clearly shows the interest of this end-of-life option. Actually, it almost offsets the whole environmental impacts of the composite manufacturing (Figure 4).

By recycling a product mainly sourced with carbon fossil fuel (Figure 1), impacts on climate change or fossil depletion can be almost completely avoided (Figure 4). For marine eutrophication, recycling allowed for a larger avoidance than the impacts of manufacturing. This is due to the use of European electricity for the injection molding of the matrix, while we use a French mix for recycling process (it impacts systematically onto this indicator).

When comparing the environmental impacts of a 1 kg composite part during its life cycle, depending on the end-oflife scenario (landfill or recycling of carbon fibers), and despite electricity consumption in the recycling process, emission of greenhouse gases may be divided by 10 (Figure 5). The environmental gain is on average about 80%, according to the ReCiPe Midpoint (H) method. For the climate change indicator, it is about 100%. This is because of the use of a French electricity country mix, which is mainly sourced by nuclear energy, which is energy that has no impacts on climate change (it does impact principally on the ionizing radiation category).

Negative impacts (for eutrophication and natural land transformation indicators) do not mean that they are "good" for the environment. This only means that it is an avoided impact; to recycle, allows for avoiding some impacts due to the manufacture stage.

Economic Validation. We recently made a market study showing that there will always be relevant uses for recycled reinforcements or for semi-products based on secondgeneration fiber, whatever their mechanical characteristics are and as long as the price remains reasonable. [START_REF] Alcimed Report | Evaluation du gisement de fibres de carbone àrecycler et de leur valorisation pour le consortium RECCO[END_REF] The integration of recycled carbon fiber is only interesting if the mechanical performance/price ratio is higher than that of glass fiber. Therefore, in light of excellent second-generation reinforcement mechanical properties, 27 this ratio should be much higher than for new carbon fibers. Thus, the feasibility of recycling will be provided if the second-generation semi-products price does not exceed 70-80% of the new ones.

■ DISCUSSION

In the present context, the use of carbon/epoxy composite is ever increasing. As indicted, these composites can be recycled by solvolysis, [START_REF] Morin; Loppinet-Serani | Near suband supercritical solvolysis of carbon fibres reinforced polymers (CFRPs) for recycling carbon fibers as a valuable resource: State of the art[END_REF] keeping good mechanical properties. [START_REF] Goedkoop | A Life Cycle Impact Assessment Method Which Comprises Harmonised Category Indicators at the Midpoint and the Endpoint Level[END_REF] Anticipating that they may soon be subjected to regulation, it is essential to show it is feasible that a composite recycling network can be set up that is both economically and environmentally favorable.

The recovery of the carbon reinforcement (which is the most environmentally impacting constituent in the composite manufacturing) by an aqueous solvolysis of the composite's matrix leads to an average gain of about 80% for all ecoindicators compared to the landfill end-of-life option.

Lastly, the remanufacturing process developed allows for obtaining a semi-product easily usable. Consequently, from an economic point of view, the mechanical performance/price ratio of the second-generation carbon fiber should be higher than that for the virgin carbon fibers or the glass reinforcement. The next step in the maturation of this technology is the development of a pilot scale facility for the recycling of carbon fibers from CFRPs using the supercritical fluid technology.

Figure 1 .

 1 Figure 1. Environmental impacts due to the carbon reinforcement (dark gray), epoxy matrix (light gray), and injection molding process (intermediate gray), while processing a 1 kg carbon/epoxy composite part. The analysis is based on Duflou et al. data.[START_REF] Duflou | Environmental impact analysis of composite use in car manufacturing[END_REF]

  Obtained by Hydrolysis in Supercritical Water. The hydrolysis of the epoxy resin matrix in supercritical water (p c = 22.1 MPa, T c = 374 °C) has been published many times as well as the alcoholysis in supercritical alcohols (methanol, p c = 8.1 MPa, T c = 239.3 °C; ethanol, p c = 6.1 MPa, T c = 240.8 °C, or still isopropanol, p c = 4.8 MPa, T c = 235.1 °C). For instance, Okajima et al. have studied the hydrolysis of epoxy resin of CFRPs in sub-and supercritical water in the temperature range between 300 and 450 °C and 25

Figure 2 .

 2 Figure 2. System boundary and life cycle stages. Dashed lines represent life cycle stages that were cut off.

Figure 3 .

 3 Figure 3. SEM image of recycled fibers after supercritical water treatment at 400 °C and 25 MPa at an ICMCB laboratory.

Figure 4 .

 4 Figure 4. Life-cycle impact assessment of the landfill of a 1 kg carbon/epoxy composite part (dark gray) compared with the reinforcement's recycling (light gray). The analysis is based on the ReCiPe Midpoint (H) method.

Figure 5 .

 5 Figure 5. Comparison of the environmental impacts of a 1 kg composite part during its life cycle, depending on the end-of-life scenario of carbon fibers (landfill in dark gray; recycling in light gray).

■ AUTHOR INFORMATION

Corresponding Author *Fax: (+33)55640006994. E-mail: g.sonnemann@ism.u-bordeaux1.fr.

Notes

The authors declare no competing financial interest.

■ ACKNOWLEDGMENTS

The authors express their thanks to the Aquitaine Regional Council for their financial support and to all the partners of the RECCO project, Innoveox, Herakles (member of Safran group), Astrium ST, and Airbus (members of AIRBUS group) and ADEME (French Environment and Energy Management Agency) for their financial support.