Camille Valdeyron

Martin Soubrier

Bruno Pereira

Arnaud Constantin

Jacques Morel

Philippe Gaudin

Bernard Combe

Anne Sophie Gremeau

Lydie Dejou-Bouillet

Jean Luc Pouly

Lydie Dejou- Bouillet

Vincent Sapin

Charlotte Oris

Florence Brugnon
email: fbrugnon@chu-clermontferrand.fr

Pr Florence

Impact of disease activity and treatments on ovarian reserve in patients with rheumatoid arthritis in the ESPOIR cohort

Keywords: Rheumatoid Arthritis, Inflammation, Fertility, Woman, Methotrexate, Ovarian

à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Rheumatoid arthritis (RA) may affect the fertility of women suffering from this pathology [START_REF] Junco | The relationship between rheumatoid arthritis and reproductive function[END_REF][START_REF] Provost | Fertility and infertility in rheumatoid arthritis[END_REF][START_REF] Ostensen | Rheumatoid arthritis: The effect of RA and medication on female fertility[END_REF][START_REF] Kay | subfertility before and after the development of rheumatoid arthritis in women[END_REF][START_REF] Katz | Childbearing decisions and family size among women with rheumatoid arthritis[END_REF][START_REF] Østensen | New insights into sexual functioning and fertility in rheumatic diseases[END_REF]]. Recent studies have shown unexplained infertility [START_REF] Eudy | Fertility and Ovarian Reserve among Women with Rheumatoid Arthritis[END_REF], with longer mean time to pregnancy [START_REF] Clowse | Effects of infertility, pregnancy loss, and patient concerns on family size of women with rheumatoid arthritis and systemic lupus erythematosus[END_REF][START_REF] Brouwer | Fertility in women with rheumatoid arthritis: influence of disease activity and medication[END_REF][START_REF] Nelson | Fecundity before disease onset in women with rheumatoid arthritis[END_REF][START_REF] Jawaheer | Time to pregnancy among women with rheumatoid arthritis[END_REF], higher prevalence of dysovulation [START_REF] Brouwer | Subfertility in rheumatoid arthritis is often unexplained or caused by anovulation[END_REF][START_REF] Flaisler | A study of ovarian function in rheumatoid arthritis[END_REF], nulliparous women [START_REF] Wallenius | Fertility in women with chronic inflammatory arthritides[END_REF], and reduced ovarian reserve [START_REF] Brouwer | Subfertility in rheumatoid arthritis is often unexplained or caused by anovulation[END_REF][START_REF] Anasti | Premature ovarian failure: an update[END_REF] compared with the general population. However, it is not clear if these women's subfertility is due to RA activity and/or their treatments. Ovarian reserve can be evaluated by measuring the serum anti-Müllerian hormone (AMH) level. AMH is a hormone produced by granulosa cells surrounding oocytes in growing follicles that have undergone recruitment in the ovary [START_REF] Durlinger | Regulation of ovarian function: the role of anti-Müllerian hormone[END_REF][START_REF] Marca | Serum antimüllerian hormone levels in women with secondary amenorrhea[END_REF][START_REF] Marca | Normal serum concentrations of anti-Müllerian hormone in women with regular menstrual cycles[END_REF]. Serum AMH level is considered to be an accurate biomarker of ovarian reserve, capable of reflecting the size of the ovarian follicular pool at reproductive age [START_REF] Tehrani | A single test of antimullerian hormone in late reproductive-aged women is a good predictor of menopause[END_REF][START_REF] De Vet | Antimüllerian hormone serum levels: a putative marker for ovarian aging[END_REF][START_REF] Sermondade | Serum antimüllerian hormone is associated with the number of oocytes matured in vitro and with primordial follicle density in candidates for fertility preservation[END_REF]. The value of this biomarker is not influenced by the period of the menstrual cycle and will reflect an early reduction of the ovarian follicle reserve [START_REF] Ferté-Delbende | Evaluation of the ovarian reserve[END_REF].

To mitigate the decrease in follicle reserve linked to gonadotoxic treatment or disease, fertility preservation is currently possible for these women through the cryopreservation of oocytes and/or ovarian tissue [START_REF] Shapira | First delivery in a leukemia survivor after transplantation of cryopreserved ovarian tissue, evaluated for leukemia cells contamination[END_REF]. When needed, the thawed oocytes may be used for assisted reproductive technology by intra-cytoplasmic sperm injection, and ovarian tissue may be grafted if infertility is confirmed. Recent studies reporting the impact of RA and its treatment on fertility were performed with a short time of follow-up [START_REF] Brouwer | Levels of serum anti-Müllerian hormone, a marker for ovarian reserve, in women with rheumatoid arthritis[END_REF] or only small patient cohort sizes [START_REF] Henes | Ovarian reserve alterations in premenopausal women with chronic inflammatory rheumatic diseases: impact of rheumatoid arthritis, Behçet's disease and spondyloarthritis on anti-Müllerian hormone levels[END_REF][START_REF] Brouwer | Levels of serum anti-Müllerian hormone, a marker for ovarian reserve, in women with rheumatoid arthritis[END_REF]. Hence, there is no consensus in the literature on indications of fertility preservation for these patients regarding RA activity and the treatments they receive.

The aim of our study was to examine the impact of RA activity and its related treatments on ovarian reserve measured by serum AMH level in a large cohort of patients for 36 months.

The information obtained in this analysis may help to better define indications of fertility preservation for these patients. and April 2005. To be included, patients had to have RA present for less than 6 months before inclusion or undifferentiated inflammatory arthritis (UA) with potential for developing RA, also for less than 6 months. UA is early arthritis that does not meet the criteria for another diagnosis, such as erythema nodosum or sacroiliitis, or that has joint characteristics and extraarticular characteristics incompatible with RA.

Patients and Methods

Patients

All the patients recruited in the ESPOIR cohort were naive to any treatment at inclusion. Only anti-inflammatory medication taken for less than 15 days before inclusion was allowed.

Patients were followed up by the same investigator every 6 months for 2 years and every year for at least 10 years. The protocol of the ESPOIR cohort was approved by the ethics committee of Montpellier, France (no. 020307). All the patients gave their signed informed consent before inclusion. One biological resources centre (Sarah Tubiana, Paris-Bichat) was in charge of centralising and managing biological data collection.

In our study, we enrolled women (n = 102) 18-37 years old who fulfilled the 2010 American College of Rheumatology/European League Against Rheumatism criteria for RA. Fifteen women without treatment who were considered to have undifferentiated arthritis (UA) with the potential to become RA and were used as a chronic inflammatory disorder control group. This group allowed us to analyse the impact of inflammation itself on ovarian reserve without confounding by any potential effect of treatment. We excluded men, and women under the age of 18 and over 37 years old. We chose the age limit of 37 years because the serum AMH level decreases significantly and physiologically at this age due to follicular atresia [START_REF] Marca | Normal serum concentrations of anti-Müllerian hormone in women with regular menstrual cycles[END_REF].We also excluded patients with ovarian surgery and premature ovarian insufficiency (POI) defined, according to the ESHRE guideline [START_REF] Webber | ESHRE Guideline: management of women with premature ovarian insufficiency[END_REF] on the management of POI 2015, by loss of ovarian activity before the age of 40 years. POI is characterized by menstrual disturbance (amenorrhea or oligomenorrhea) with raised gonadotropins and low estradiol.

Measurements

The primary outcome was the global evolution of the serum AMH level (ng/mL), adjusted for age, over 36 months of follow-up in patients with RA. The measurement of AMH level was performed using the electrochemiluminescence method with a Cobas® e411 analyser (Roche Diagnostics, Meylan, France). The threshold of <1 ng/mL was considered to be a proof of severe decrease in ovarian reserve. The serum AMH assay was performed on patients from the ESPOIR cohort at baseline (T0), 6, 12, 24, and 36 months after inclusion (T6, T12, T24, and T36, respectively). All the results were expressed as the median value of AMH level (ng/mL). Blood samples were available for 117 patients at T0, 111 at T6, 106 at T12, 99 at T24, and 82 at T36 months of follow-up.

We analysed the impact of RA activity and inflammation on the serum AMH level. Disease activity was scored using the Disease Activity Score in 28 joints (DAS28 CRP) and biological parameters of inflammation: C-reactive protein in mg/L (CRP) and erythrocyte sedimentation rate in mm/h (ESR). Anti-CCP antibodies were measured from thawed frozen sera using Immunoscan CCplus®, Eurodiagnostica, Arnhem, The Nederlands.

Finally, we evaluated the potential impact of treatments used on AMH level. The treatments used were methotrexate (MTX) only or with other treatments. The focus on MTX was justified by its high frequency of use because it is the first line of treatment used in RA according to the EULAR and SFR recommendations (29) [START_REF] Smolen | EULAR recommendations for the management of rheumatoid arthritis with synthetic and biological disease-modifying antirheumatic drugs: 2016 update[END_REF]. To study the potential effects of the treatments, we analysed four groups of patients according to their treatment (Fig. 2):

MTX only (n = 27), MTX and other treatments (n = 56) (hydroxychloroquine, leflunomide, sulfasalazine, etanercept, adalimumab, rituximab, and infliximab), treatment without MTX (n = 19) (leflunomide, sulfasalazine, and hydroxychloroquine), and a control group without treatment (n = 15). We also studied the impact of corticosteroid on AMH level, with 69 exposed patients and 37 patients without corticosteroid treatment.

Statistical analysis

Statistical analyses were performed using Stata software (version 13, StataCorp, College Station, US). All tests were two-sided, with a Type I error set at 0.05. Continuous data were expressed as median (interquartile range) according to the statistical distribution. The assumption of normality was assessed using the Shapiro-Wilk test.

Concerning baseline data, the comparisons were conducted using (i) chi-square and Fisher's exact tests, as appropriate, for categorical data, and (ii) analysis of variance (ANOVA), or the Kruskal-Wallis test when the assumptions of ANOVA were not met, for continuous parameters. When appropriate (omnibus P-value less than 0.05), a post-hoc test for multiple two-by-two comparisons was performed: Tukey-Kramer after ANOVA and Dunn post Kruskal-Wallis. To study the relationships between continuous parameters, Pearson or Spearman correlation coefficients were estimated. Random-effects models (more precisely, linear mixed models) were performed to compare the evolution of the AMH level between treatment groups, (i) studying the following fixed effects: group, time-points, and their interaction, and (ii) taking into account the between-and within-patient variability (considering the subjects as random-effect). Sidak's type I error correction was applied to perform multiple comparisons. The normality of residuals was analysed using the Shapiro-Wilk test and a logarithmic transformation of the AMH level was then conducted to assess the normality. Multivariable analyses were then carried out with adjustment for age and DAS28 at inclusion, as the treatment groups differed in these parameters at inclusion. A sensitivity analysis was performed, considering the AMH level as a normal/non-normal dichotomous parameter. The age-related threshold reported in the literature was applied. To compare the evolution of the AMH level between treatment groups, a generalized linear mixed model with a logit link function was used.

Results

Baseline demographics and disease activity

Study population details are shown in Table 1. In all, 117 women from the ESPOIR cohort were included. They were 29.9 (24.8-33.5) years old. Among these enrolled women, 102 (90%) fulfilled the 2010 ACR/EULAR RA classification criteria and were diagnosed as RA patients. They began their treatment on the day of inclusion (T0). Fifteen (10%) women, diagnosed as having UA without treatment, were included in the chronic inflammatory disorder control group. These women were 28.2 (25.0-33.8) years old.

Median body mass index (BMI) was 22.2 (20.5-25.0) kg/m 2 . Sixty-seven patients (51%) were smokers. As shown in Table 1, there were no differences in these characteristics between the evaluated patients and the control group.

Rheumatoid factor (RF) was positive in 43 (33%) patients and anti-citrullinated protein antibodies (ACPA) were positive in 45 (35%) patients. There was a significant difference in both these parameters between the subgroup of patients treated with MTX and the other groups without MTX (p<0.05). There was no difference in serum AMH level between patient with ACPA-positive patients and ACPA-negative patients, with respective values of 1.52 (0.63-3.26) ng/mL and 1.72 (0.97-2.81) ng/mL (P=0.64). The median serum CRP level was 10 (5-24) mg/L and the median ESR level was 20 (10-35) mm/h. There was no significant difference in CRP and ESR levels at T0 between any groups of patients.

The patients' median value of DAS28 was 5.0 (4.3-5.8), corresponding to moderate disease activity. The DAS28 score of the subgroup of patients treated with MTX was higher at 5.5(4.5-6.2) corresponding to high RA activity compared with the other groups (P = 0.001), and a higher level of inflammation and disease activity for these patients. The Health Assessment Questionnaire (HAQ) score was significantly higher in the MTX group (1.25[0.88-1.63]) compared with the other groups (P < 0.02).

Impact of disease activity and inflammation on serum AMH level

The analysis of serum AMH level over the time and as function of age describes a gradual decrease (Fig. 1). There was no significant difference between AMH level at T12, T24, and T36 compared with T0 (P > 0.05). We observed a lower AMH level at T0 compared with T6

(1.70 ng/mL vs. 2.20 ng/mL; P = 0.017) in the group of 102 patients with treatment. Control group serum AMH levels measured at the different times during the 36-month follow-up period did not differ significantly.

At T0, the serum AMH level was negatively correlated with the DAS28 score (r = -0.12; P = 0.15) and with the biological parameters of inflammation (CRP, r = -0.26; ESR, r = -0.19; P = 0.002) (Fig. 3).

There was a statistically significant decrease in RA activity (measured by DAS28 score; r = -0.28; P = 0.003) and a trend for inflammation (measured by serum CRP level; r = -0.17; P = 0.06) which was correlated with the increase in AMH level (1.70 ng/mL vs. 2.20 ng/mL) between T0 and T6 for all the patients evaluated, whatever treatment they received. The control of inflammation was correlated with an improvement in AMH level. These variations were not observed in the control group (r = 0, P = 0.99).

There was no statistical difference between groups regarding the evolution of AMH level between T6, T12 (P = 0.09), T24 (P = 0.35), and T36 (P = 0.61).

Impact of treatments on serum AMH level during 36 months of follow-up

At T0, we found a lower serum AMH level in the "MTX only" patient group (n = 27) compared with the "without MTX" patient group (n = 19; P = 0.002) and the control group (n = 15; P > 0.05). The DAS28 score was higher in the "with MTX" patient group (DAS28 at 5.5 (4.5-6.2)) compared with the "without MTX" patient group (DAS28 at 4.6 [3.8-5.5 ]), P=0.001.

Although the mean AMH level tended to decrease during the follow-up, we did not observe a significant difference between the serum AMH levels at T12, T24, and T36 compared with T0 (P > 0.05), whatever the treatments (MTX or not) in the patients evaluated, even after adjustment for the DAS28 score. Indeed, the serum AMH level evolved from 1.13 ng/mL to 1.14 ng/mL at T36 in the MTX-only group (n = 27), from 1.52 to 1.34 ng/mL at T36 in the MTX with or without other treatments (n = 83) group, and from 1.72 ng/mL to 1.62 ng/mL at T36 in the without MTX group (n = 19), versus from 1.62 ng/mL to 0.94 ng/mL in the control group (n =15) (Fig. 2).

There was no difference in the evolution of the serum AMH level (Fig. 1) over time between RA patients under treatment (n = 102) and the control group (n = 15; P = 0.13) and between groups of patients who fulfilled (n=102) or did not fulfil (n=15) the ACR/EULAR 2010 criteria (P>0.05).

In regard to corticosteroid treatments, we did not find a difference in serum AMH levels between users (n = 69) and non-users (n = 37) at any time of the comparison, even with doses higher than 7.5 mg per day (n = 6).

Discussion

By following this cohort of women of reproductive age with early RA for a period of 36 months, we were able to clearly show that RA activity was correlated with a decrease in serum AMH level. However, AMH level as a function of age in these patients is similar to the physiological concentration of AMH over time reported by Kelsey et al. (31) in women without RA. In addition, we did not observe any impact of treatments used to treat RA on AMH level of these women.

This study enabled us to show clearly to show the impact of RA activity on serum AMH level early in the disease. This decline in serum AMH level is linked to the decrease in ovarian reserve we observed when diagnosing RA before starting treatment. The strength of our study lies in the size of the ESPOIR cohort and the long follow-up period of 36 months. No previous study has carried out as lengthy a longitudinal follow-up of an equally large number of patients. We were only able to evaluate ovarian reserve on the basis of serum AMH level because ovarian echography had not been performed for the patients included in the study.

The impact of RA on female fertility remains controversial. Henes et al. [START_REF] Henes | Ovarian reserve alterations in premenopausal women with chronic inflammatory rheumatic diseases: impact of rheumatoid arthritis, Behçet's disease and spondyloarthritis on anti-Müllerian hormone levels[END_REF] showed that patients with RA had a lower AMH level, a lower pregnancy rate (Del Junco et al. [START_REF] Junco | The relationship between rheumatoid arthritis and reproductive function[END_REF], and a longer waiting time to pregnancy (8,9) compared with the general population. Other authors focused more particularly on the effect of MTX, since it is a first-line treatment for RA. Two studies pointed to a reduction of ovarian reserve after exposure to MTX [START_REF] Kelsey | A validated model of serum anti-müllerian hormone from conception to menopause[END_REF][START_REF] Jancin | Methotrexate may curb girls' ovarian function[END_REF]. Nonetheless, these observations were not confirmed by Brouwer et al. [START_REF] Brouwer | Levels of serum anti-Müllerian hormone, a marker for ovarian reserve, in women with rheumatoid arthritis[END_REF] in a follow-up of RA patients over a 6-months period of exposure to MTX. A similar result was observed in a recent study of Eudy et al. in 2019 [START_REF] Eudy | Fertility and Ovarian Reserve among Women with Rheumatoid Arthritis[END_REF], with no impact of MTX on AMH level. In our study, we did not find an impact of MTX during the 36-month follow-up period in the patients exposed only to MTX (n = 27) or in those exposed to MTX in addition to other treatments (n = 56). The diversity of groups of treatments in our study is explained by the fact that the cohort of patients is multicentred and by the difficulty of clearly determining the diagnosis of RA in the framework of early arthritis. Such diagnosis involves a variety of recommendations and diverse clinical rheumatologic practices, as shown by Benhamou et al. (34). This likely explains our study's possible lack of power in its evaluation of the effects of treatments, due to the lower number of patients in each group after adjustment. However, it is noteworthy that a large proportion of the treatments used are not known to affect fertility-this reduces the risk of interference in the interpretation of our results with, in particular, leflunomide, hydroxychloroquine, and sulfasalazine. Only a severe impact would justify preserving fertility. Our study therefore confirms results published previously that do not recommend preserving fertility for patients undergoing treatment with immunosuppressive doses of MTX [START_REF] Choux | Immunosuppressive therapy and fertility preservation: Indications and methods[END_REF]. However, it is important to note that since the formation of this cohort new and effective RA treatments have been developed. It will be interesting to analyse the impact of new generation of treatments on AMH in a further study.

Regarding corticosteroid therapy, no significant difference in serum AMH level was observed between the exposed and unexposed patients. This observation does not confirm the data in the literature. Brouwer et al. [START_REF] Brouwer | Fertility in women with rheumatoid arthritis: influence of disease activity and medication[END_REF] highlighted a link between corticosteroid therapy and a reduction in fertility in 66% of the 245 women in the Pregnancy-Induced Amelioration of RA (PARA) cohort composed of patients following corticosteroid therapy at a dose of 7.5 mg/day. When patients followed corticosteroid therapy at a dose higher than this daily dose, an increase in the waiting time to pregnancy was observed in comparison with a control group. In our study, the number of patients exposed to this dose was very low, which perhaps explains the observed absence of difference.

The analysis of AMH levels according to ACPA status did not show any significant difference contrary to the previous study in 2019 by Brouwer et al. (36) which demonstrated lower AMH levels in ACPA-positive patients compared with healthy women as a control group. In our study the control group consisted of patients with UA. Moreover, we applied the same assay to measure serum AMH level in patients and control groups. This may explain this lack of difference in our study.

No correlation was found in previous studies between serum inflammatory parameters, disease activity, and AMH level. Brouwer et al. [START_REF] Brouwer | Reduced Ovarian Function in Female Rheumatoid Arthritis Patiens trying to conceive[END_REF] did not show any association between disease activity and serum AMH level in a recent study. However, their study was based on a population whose disease activity was under control and who had participated in a pregnancy project. Our study clearly shows a correlation between RA activity and AMH level. At inclusion, inflammation linked to RA was not controlled by treatment, with RA activity evaluated as moderate or high as a function of group, according to the DAS28 score. The organisation of treatment from the time of diagnosis of RA makes it possible to see its probable beneficial impact on disease evolution at 6 months, with a significant improvement in serum AMH level whatever the treatment group. An increased prevalence of infertility of unknown aetiology was shown by Brouwer et al. [START_REF] Brouwer | Subfertility in rheumatoid arthritis is often unexplained or caused by anovulation[END_REF] in these patients. Our data appear to show that it might be an effect of inflammation on the ovarian reserve.

RA activity, which reflects the magnitude of articular inflammation, generally appears to have a significant impact on serum AMH level. Previously, it was only possible to draw a link between increased waiting time to pregnancy and the inflammatory activity of RA [START_REF] Ostensen | Rheumatoid arthritis: The effect of RA and medication on female fertility[END_REF]. The hypothesis used to explain this association was potential sexuality disorders linked to a fall in libido related to inflammation-induced pain. Pikwer et al. showed an association between the occurrence of the menopause before the age of 45 and the development of RA, with an odds ratio (OR) of 2.42, 95% CI 1.32 to 4.45 [START_REF] Bergström | Early menopause is an independent predictor of rheumatoid arthritis[END_REF]. The authors hypothesized a link between the inflammation present before the diagnosis and the reduction of ovarian reserve, leading to premature menopause. The same observation was made by Del Junco et al., who demonstrated that female fertility was reduced before the diagnosis of RA was made [START_REF] Junco | The relationship between rheumatoid arthritis and reproductive function[END_REF]. In our study also, it was at T0 that we found the lowest AMH level. In comparison with other models of chronic systemic diseases such as Crohn's disease, it has been shown that serum AMH values are also clearly reduced in patients with intense disease activity compared with a control group of patients in remission [START_REF] Brouwer | Reduced Ovarian Function in Female Rheumatoid Arthritis Patiens trying to conceive[END_REF]. Inflammation therefore appears to play a crucial role in fertility in RA.

Indeed, proinflammatory cytokines such as interleukin (IL) and TNF can harm the vascular system, possibly explaining the link between the inflammatory activity of a disease and damage to the ovarian reserve [START_REF] Choux | Immunosuppressive therapy and fertility preservation: Indications and methods[END_REF][START_REF] Bergström | Early menopause is an independent predictor of rheumatoid arthritis[END_REF]. Paradisi et al. [START_REF] Paradisi | High cytokine expression and reduced ovarian reserve in patients with Hodgkin lymphoma or non-Hodgkin lymphoma[END_REF] studied the potential link between cytokines and the ovarian reserve in patients with Hodgkin's and non-Hogkin's lymphoma before the administration of chemotherapy. A negative correlation was found between AMH level and IL6 in these patients [START_REF] Paradisi | High cytokine expression and reduced ovarian reserve in patients with Hodgkin lymphoma or non-Hodgkin lymphoma[END_REF]. In an animal model, IL-1 was implicated in the regulation of the ovarian reserve through an increase of serum AMH level, with a better response of the ovaries to gonadotropins in IL-1-beta-deficient mice [START_REF] Shaish | Interleukin-1 deficiency prolongs ovarian lifespan in mice[END_REF]. It would therefore be interesting to extend our study by measuring cytokine levels in the serum and/or follicular fluid to analyse their potential influence and better understand how inflammation affects the ovarian reserve of patients with RA. One might imagine that reducing inflammation would remove this "shock" to the production of AMH by the granulosa cells. Accordingly, it would be interesting to know whether the granulosa cells can be destroyed by long exposure to inflammation (driven by cytokines), since that would clearly link inflammation to a reduced ovarian reserve.

In conclusion, our study highlights the harmful impact of RA activity on the serum AMH level, an indicator of the ovarian reserve. We did not demonstrate the impact of treatments.

Our study therefore tends to show that patients with RA do not need to preserve their fertility.

Nevertheless, it appears important to thoroughly control the inflammatory activity of RA from the time it is considered as a diagnosis, in order to limit the premature impact of inflammation on the ovarian reserve. In the future, it would be interesting to measure interleukin levels in follicular fluid to better understand the physiopathology of the damage incurred. The relationship between age of the patients and serum AMH level was analysed using fractional-polynomial prediction (blue line). Red line corresponds to the 95% confidence interval. The serum AMH level at the diagnosis of the disease is negatively correlated with the DAS 28 score (r= -0.12; p=0.15).

  Patient data and blood sample analyses were obtained from the French prospective multicentred observational cohort ESPOIR (27) (Etude et Suivi des Polyarthrites Indifférenciées Récentes) sponsored by the French Society of Rheumatology (FSR), in December 2002. The ESPOIR cohort is a French database for various investigations (diagnostic, prognostic, economics, etc.). The study included 813 patients aged 18-70 years with early arthritis and attending 14 rheumatology centres in France between December 2002

Fig. 1

 1 Fig. 1 Evolution of serum AMH level during 36 months of follow-up according to the age of

Fig. 2

 2 Fig. 2 Evolution of serum AMH level for each treatment

Fig. 3

 3 Fig. 3 Correlation between serum AMH level and parameter of inflammation :

Acknowledgements:

An unrestricted grant from Merck Sharp and Dohme was allocated for the first 5 years. Two additional grants from INSERM were obtained to support part of the biological database. The French Society of Rheumatology, Pfizer, Abbvie, and Lilly also supported the ESPOIR cohort study. We also wish to thank Nathalie Rincheval (CHU Montpellier and EA 2415) who did expert monitoring and data management, and all the investigators who recruited and followed the patients (F. Berenbaum, Paris-Saint Antoine; M.C. Boissier, Paris-Bobigny; A. Cantagrel, Toulouse; B. Combe, Montpellier; M. Dougados, Paris-Cochin; P. Fardellone and P. Boumier, Amiens; B. Fautrel, Paris-La Pitié; R.M. Flipo, Lille; Ph. Goupille, Tours; F. Liote, Paris-Lariboisière; O. Vittecoq, Rouen; X. Mariette, Paris Bicetre; P. Dieude, Paris Bichat; A. Saraux, Brest; T. Schaeverbeke, Bordeaux; and J. Sibilia, Strasbourg). We finally wish to thank A. Vega and C.Lambert (Biostatistics Unit (DRCI)), CHU Clermont-Ferrand for our contribution. This work was supported by Roche Diagnostics with the generous gift of the AMH assays.

Conflicts of interest

No conflicts of interest.


Jean-Philippe Argaud
email: jean-philippe.argaud@edf.fr

Sibo Cheng
email: sibo.cheng@edf.fr

Bertrand Iooss
email: bertrand.iooss@edf.fr

Didier Lucor
email: didier.lucor@limsi.fr

Angélique Ponçot
email: angelique.poncot@edf.fr

Iterative methods for improving error covariance modeling in variational assimilation

come

ABSTRACT

Variational data assimilation methods are well known and widely applied in geophysical domains to problems affected by uncertainty of short-and long-term model predictions. The strength of this approach is to fuse available information in order to find a compromise between background model predictions and observations where the associated weights are provided by error covariance matrices. A remarkable difficulty in data assimilation is the lack of information for background error covariance modeling. In many industrial applications, due to the absence of historical observations/predictions, the background matrix modelling remains empirical relying on some form of expertise and imposed physical constraints. This can be problematic especially for short term prediction or static reconstruction when the ensemble methods become inappropriate.

Great effort has been made to diagnose the covariance matrices modelling a posteriori. For instance, the meteorology community (e.g. Météo-France) has been a strong contributor to this topic. Several algorithms and analysis were proposed, especially the sucessful iterative method introduced in [START_REF] Desroziers | Diagnosis and adaptive tuning of observation-error parameters in a variational assimilation[END_REF] which adjusts sequentially the error covariance ratios in multivariate systems. Recent efforts are also carried out to apply Desroziers iterative methods in local subspaces, which could make this method more flexible. In this work, under the assumption of a good knowledge of the observation error covariance, we have developed two sequential adaptive methods for improving the assimilation result and rebuilding output error covariance matrices. These two methods formalize existing practice, consisting in repeating several times the assimilation procedure with the same observations. In the framework of twin experiments with both state-independent and state-dependent noises, we show that these methods, as well as the method of Desroziers are helpful in terms of decreasing the error expectation and the innovation quantity (observationanalysis gap).

We have then assessed these methods in a context of a data assimilation chain developed for the purpose of a hydrological problem whose objective is to improve the prediction of the outflow of a watershed.


Siham El Garroussi
email: garroussi@cerfacs.fr

Matthias De Lozzo

Sophie Ricci
email: ricci@cerfacs.fr

Didier Lucor
email: didier.lucor@limsi.fr

Nicole Goutal
email: nicole.goutal@edf.fr

Cédric Goeury
email: cedric.goeury@edf.fr

Sébastien Boyaval
email: sebastien.boyaval@enpc.fr

Siham El Garroussi

Matthias De Lozzo
email: matthias.delozzo@irt-saintexupery.com

Uncertainty quantification in a two-dimensional river hydraulic model

Keywords: Open-channel flow, Sensitivity analysis, Surrogate model, Gaussian process, Monte Carlo method

the spatially distributed water level with respect to uncertain stationary flow to the model and friction coefficients. The quality of the surrogate was assessed on a validation set, with small root mean square error and a predictive coefficient equal to 1. Sobol' sensitivity indices are computed and enhance the high impact of the input discharge on the water level variation.

INTRODUCTION

Flood inundation models are central components in any flood risk analysis system as they transform the bulk discharge outputs from flood-frequency analyses or rainfall-runoff models into distributed predictions of flood hazard in terms of water depth, inundation extent and flow velocity. Predictions may be dynamic in time and can be derived from a range of codes which vary in complexity from non-model approaches, such as fitting a planar surface to digital elevation data, through to numerical solutions of fluid dynamics equations derived from considerations of mass and momentum conservation. Whilst such models are parsimonious in terms of their data requirements and number of unconstrained parameters relative to other environmental physics, their underlying equations may be non-linear. Moreover, the data sets that they do require may be subject to complex, but poorly known errors that may vary markedly in time and space. As a consequence, considerable research has, in recent years, sought to understand and better estimate these uncertainties in order to improve flood risk analysis. Typically, uncertainties in hydrodynamics models stand are classified as: parametric (or epistemic) uncertainty, arising from incomplete knowledge of the correct settings of the models parameters; input data uncertainty, arising from incomplete knowledge of the true value of the initial state and forcing, usually linked to the aleatory nature of the physics; and structural uncertainty, which is the failure of the model to represent the system, even if the correct parameters and inputs are known. Together, these three components represent a complete probabilistic description of the informativeness of the model for the underlying system. But in practice, all are extremely challenging to specify.

In this study, we consider both epistemic and aleatory uncertainties by investigating the effect of two uncertainty sources on water level calculation for extreme flood event, respectively the roughness coefficient and the upstream discharge. On the one hand, the estimation of the roughness is difficult because it is a lumped parameter that mostly reflects the flow resistance of the river. Since the roughness coefficient has an extensive effect on flow analysis of a river, including computation of the water level and velocity, its accurate estimation is important for prediction of the water level during flooding. Because of its importance, various efforts have been made to quantify the roughness coefficients of rivers in an objective manner. Among them, an element-based method [START_REF] Cowan | [END_REF] and empirical equations that relate the roughness coefficient either to bed material [START_REF] Strickler | Beitrge zur Frage der Geschwindigkeitsformel und der Rauhigkeitszahlen fur Strme, Kanle und Geschlossene Leitungen[END_REF] or to relative depth [START_REF] Charlton | Hydraulics Research Station. The hydraulic geometry of some gravel rivers in Britain[END_REF] are representative. However, owing to the diversity and irregularity of natural rivers, prediction of the roughness coefficient for a specific river reach using these methods is not simple. Thus, until now, field measurements have been made either to directly estimate the roughness coefficient [START_REF] Coon | Estimates of roughness coefficients for selected natural stream channels with vegetated banks in New York[END_REF] or to provide references [2,[START_REF] Hicks | Roughness characteristics of New Zealand rivers[END_REF]. However, there remain uncertainties whether using the methods referred to above or using field measurements. From a practical viewpoint, water level and discharge as variables computed by numerical modeling are influenced by uncertainty in estimating the roughness coefficient. Conducting simulation of dam breakage flow for the Teton Dam, [START_REF] Fread | BREACH: an Erosion Model for Earthen Dam Failures (Model description and User Manual)[END_REF] showed that variation in calculated flood flow water depth was less than 5% with a 20% change in the roughness coefficient. He therefore argued that even if uncertainty in Stricklers roughness coefficient is large, its effect on the water depth might be reduced considerably in the process of computation. These conclusions should be deeply investigated in the context of flood simulation, on our own river test case, characterized by long homogeneous friction zones calibrated in high flow. On the other hand, inundation models require the specification of boundary conditions, which are typically the greatest source of aleatory and epistemic uncertainty when simulating the annual exceeding probability of inundation. Flow at the upstream boundary of the river is often the most important boundary condition, although most applications will require (unless using a kinematic solution) or benefit from downstream-level boundaries (e.g. tidal reaches). In locations where they are available, gauging stations are typically the most accurate source of river flow and level data. However, the ratings at these stations, that convert observed levels to flows, are usually based on low to medium flow observations, necessitating an uncertain extrapolation of the rating for high flows. Rating errors may be especially large when flow is out of bank. Where gauging stations are not available or spatially sparse, hydrological models can be used to simulate upstream discharges. However, despite much effort, rainfall-runoff models are still very uncertain, especially where calibration/validation data are lacking.

Subsequently, once the sources of uncertainties have been identified, they must be propagated in the model. The Monte Carlo (MC) methods are the most common techniques used for uncertainty propagation (UQ) [START_REF] Ghanem | Springer. Handbook of Uncertainty Quantification[END_REF]. This framework allows to estimate standard statistics on the model output, e.g. expectation, standard deviation, quantiles or probabilities of exceeding a given threshold. It also makes it possible to estimate sensitivity indices representing the shares of output uncertainty attributable to the different uncertain input parameters, e.g. Sobol' indices where output uncertainty is measured in terms of variance [START_REF] Sobol | [END_REF]. MC is simple and highly adapted to massively parallel computational resources. Yet, its convergence is slow as it scales inversely to the square root of the sample size and its cost gets prohibitive for expensive models. In this respect, surrogate models such as Gaussian process model (GP) have received tremendous attention in the last few years, as it allows one to replace the original expensive model by a surrogate which is built from an experimental design of limited size [START_REF] Rasmussen | Gaussian Processes for Machine Learning[END_REF]. Then the surrogate can be used to compute the UQ study in negligible time. In particular, [START_REF] Roy | Stochastic Environmental Research and Risk Assessment. Comparison of polynomial chaos and Gaussian process surrogates for uncertainty quantification and correlation estimation of spatially distributed open-channel steady flows[END_REF] have shown that, for a 1D hydraulic model, on the Garonne river section that we consider and stationary flow, it features similar performance to estimate statistics by Monte-Carlo random sampling when friction and input forcing uncertainties are taken into account. The accuracy of the water level correlation matrix and sensitivity Sobol' indices estimated with the GP surrogate was assessed with respect to a classical MC estimate based on a large data set. This article is a reference for us because it involves the same section of Garonne river, the same types of uncertain variables (friction and upstream flow) and the same family of surrogate model as those considered in our work. Our work can be seen as an extension to two-dimensional flow modelling and floodplain characterization.

The present study extends the surrogate model approach in hydraulics to 2D modeling taking into account the dynamics of the flood plain. Section 2 presents the hydrodynamics solver TELEMAC-2D, the Garonne test case used in this article and the associated uncertainties. Section 3 presents the GP surrogate strategy based on the reduction of the dimension of the output space with a Proper Orthogonal Decomposition (POD). This section also presents the metrics used to assess the quality of the surrogate along with the sensitivity indices based on output variance decomposition. Results are presented in Section 4. Conclusions and perspectives are finally given in Section 5.

MODEL: TWO-DIMENSIONAL FLOW OF THE GARONNE RIVER UNDER UN-CERTAINTY 2.1 Physical model

The Shallow Water Equations (SWE), also called depth-averaged free surface flow equations, are commonly used in environmental hydrodynamics modelling [START_REF] De Saint-Venant | Théorie du mouvement non-permanent des eaux, avec application aux crues des rivières et à l'introduction des marées dans leur lit[END_REF]. They are derived from the Navier-Stokes equations [START_REF] Sohr | The Navier-Stokes Equations[END_REF] and express mass and momentum conservation averaged in the vertical dimension. The non-conservative form of the equations are written in terms of the water depth (h) and the horizontal components of velocity (u and v):

∂h ∂t + div (hu) = 0 (1)
∂u ∂t + u ∂u ∂x + v ∂u ∂y = -g ∂H ∂x + F x + 1 h div hν e --→ grad (u) (2)
∂v ∂t + u ∂v ∂x + v ∂v ∂y = -g ∂H ∂y + F y + 1 h div hν e --→ grad (v) (3)
where:

       F x = - g K s 2 u √ u 2 + v 2 h 4 3 - 1 ρ w ∂P atm ∂x + 1 h ρ air ρ w C DZ U w,x U 2 w,x + U 2 w,y F y = - g K s 2 v √ u 2 + v 2 h 4 3 - 1 ρ w ∂P atm ∂y + 1 h ρ air ρ w C DZ U w,y U 2 w,x + U 2
(h = H -z f if z f [m]
is the bottom level), u and v [m.s -1 ] are the horizontal components of velocity and ν e [m 2 .s -1 ] is the water diffusion coefficient. div and --→ grad are respectively the divergence and gradient operators.

To solve the system of equations (Eq. ( 1) to Eq. ( 3)), initial conditions h(x, y, t = 0) = h 0 (x, y); u(x, y, t = 0) = u 0 (x, y); v(x, y, t = 0) = v 0 (x, y) are provided along with boundary conditions (BC) at surface, at bottom and at upstream and downstream frontiers h(x BC , y BC , t) = h BC (t).

Study area

The study area extends over a 50 km reach of the Garonne river (France) between Tonneins (upstream), downstream of the confluence with the river Lot, and La Réole (downstream) (see Figure 1). This part of the valley was equipped in the 19 th century with infrastructure to protect the Garonne flood plain from flooding such as that occurred in 1875. A system of longitudinal dykes and weirs was progressively constructed after that flood event to protect the floodplains, organize submersion and flood retention areas. Protections on the Garonne river form a system of successive storage areas for the flood plain beyond the dikes.

Uncertainty characterization

The hydraulic variables are discretized on an unstructured triangular mesh over the twodimensional study area. We note h the vector of the water level over the p = 41416 nodes of • The roughness coefficient defined according to 4 areas. Indeed, T2D was calibrated for high flow in [START_REF] Besnard | Comparaison de modèles 1D à casiers et 2D pour la modélisation hydraulique d´une plaine d´inondation -Cas de la Garonne entre Tonneins et La Réole[END_REF] using steady-state water surface profiles at high discharge, from bankfull discharge in the main channel (2 500 m 3 .s -1 ) to bank-full discharge in the overbank flow channel between dykes. In the floodplains, the roughness coefficient K s,1 is selected as an area with cultivated fields all around the river with a Strickler coefficient of 17 m 1/3 .s -1 . Classically, according to the available expert knowledge, the friction coefficient is contained in an interval bounded by physical values depending on the roughness of soil material.

For the main channel, the Strickler roughness coefficient was split into three different areas:

from Tonneins to upstream of Mas d'Agenais, K s,2 : 45

m 1/3 .s -1 , -from upstream of Mas d'Agenais to upstream of Marmande, K s,3 : 38 m 1/3 .s -1 ,
from upstream of Marmande to La Réole, K s,4 : 40 m 1/3 .s -1 .

The distribution of Strickler roughness coefficient is chosen uniform and the interval is set to cover the range of calibration values.

• The upstream discharge is assumed to follow a Gaussian distribution centered around the thousand return period 8 490 m 3 .s -1 with a standard deviation of 700 m 3 .s -1 . The study is thus focused on extreme flood events that activate the flood plains.

Tab. 1 summarizes the considered input uncertainties.

Variable Distribution Units Q U[2 500, 10 000] m 3 .s -1 K s,1 U[5, 20] m 1/3 .s -1 K s,2 U[40, 50] m 1/3 .s -1 K s,3 U[33, 43] m 1/3 .s -1 K s,4 U[35, 45] m 1/3 .s -1
Table 1: Distributions of the input variable uncertainties.

Computing environment

In this work, hydrodynamic is provided using TELEMAC-2D (T2D) depth-averaged hydrodynamic model 1 . It solves the SWE in two dimensions with an explicit first-order time integration scheme, a finite element scheme and an iterative conjugate gradient method. In each point of the mesh, T2D gives the water depth and the vertically average horizontal velocity field [START_REF] Hervouet | Hydrodynamics of Free Surface Flows: Modelling with the finite element method[END_REF]. T2D was developed initially by the National Hydraulics and Environment Laboratory (LNHE) of the Research and Development Directorate of EDF, and is now managed by a consortium. The software comes with an API to modify the values of the uncertain parameters in a nonintrusive way.

The surrogate model construction as well as the sensitivity analysis was carried out using the BATMAN-Open-TURNS (BATMAN-OT) library2 . This library (developed at CERFACS and CECILL-B licensed) provides a convenient, modular and efficient framework for design of experiments, surrogate model and uncertainty quantification [START_REF] Roy | BATMAN: Statistical analysis for expensive computer codes made easy[END_REF]. It relies on open source python packages dedicated to statistics (openTURNS3 [START_REF] Baudin | Handbook of Uncertainty Quantification[END_REF] and scikit-learn4 [23]). It also implements advanced methods for resampling, robust optimization and uncertainty visualization.

In terms of infrastructures, CERFACS's cluster, Nemo, has been used to run T2D simulations. The Nemo cluster includes 6,912 cores distributed in 288 compute nodes. The ECU power peak is 277 Tflop/s. On this architecture, simulating the river and flood plain dynamics for the test case presented in Sect. 2.2 over 3 days, takes about 6 minutes on 15 cores.

AN EFFICIENT UQ FRAMEWORK FOR COSTLY TWO-DIMENSIONAL SIMU-LATOR

The Monte Carlo (MC) framework is the most common framework used for uncertainty quantification, due to its simplicity and good statistical results. It is theoretically applicable whatever the complexity of the deterministic model or the desired statistical estimator. However, the required sample size increases squarely with the estimator accuracy and makes this approach rather impracticable when the computational cost of each run of the model, like T2D, is non negligible. One way to lower the computationally demanding is to replace, on one side, the T2D model by a surrogate model [START_REF] Lozzo | [END_REF], on the other side, the pure random sampling by alternative sampling methods such as the Latin Hypercube sampling approach [START_REF] Goeury | Uncertainty quantification on a real case with Telamac-2D[END_REF].

Build a surrogate model with spatial output

The surrogate model based on Gaussian process regression (GP) [START_REF] Rasmussen | Gaussian Processes for Machine Learning[END_REF] has been adopted in the following. We have chosen this metamodel to the detriment of others for two reasons. The first one is its small number of hyperparameters: about one per input parameter. The second one is that it provides a measurement of its model error which would be of interest in our future work, for optimization and data assimilation problems, based on methods as expected improvement [START_REF] Jones | [END_REF]. While surrogate models offer a low cost alternative to costly models, their formulation is challenging in high dimension for inputs and outputs. In the present case, the size of the uncertain input space is small and resumes to 5 scalars. Yet, the quantity of interest is 2D and discretized over more than 41 000 points. The output space is reduced in order to limit the cost of the surrogate formulation and the spatial coherence of the later, using a Proper Orthogonal Decomposition (POD) strategy [START_REF] Ramsay | Functional data analysis[END_REF]. POD is a post-processing technique that takes a given set of data and extracts basis functions, that contain as much "energy" as possible. The meaning of "energy" depends on which kind of POD is used [START_REF] Cordier | Proper Orthogonal Decomposition: an overview[END_REF]. Here, only POD based on snapshot method [START_REF] Siade | Water resources research. Snapshot selection for groundwater model reduction using proper orthogonal decomposition[END_REF] is considered.

We propose to build a surrogate model combining POD and GP surrogate model. So we call it "POD+GP surrogate model". The corresponding algorithm is presented as follows:

1. build a learning dataset D l = x (i) , h (i) 1≤i≤N l of size N l where the design of experiments

x (i) 1≤i≤N l
is a Latin hypercube sample (LHS) [START_REF] Mckay | Technometrics. A comparison of three methods for selecting values of input variables in the analysis of output from a computer code[END_REF] with x = (Q, K s,1 , K s,2 , K s,3 , K s,4 ) and h (i) is the water level computed by T2D over the mesh at x (i) ; 2. decompose the sampled output vector h by achieving a POD on the centered output learn-

ing matrix H = h (i) j -N -1 l N l k=1 h (k) j 1≤i≤N l 1≤j≤p
and derive the most significant components; then, any sampled local water level h (i) can be approximated by a weighted sum of these components where weights depend on the input vector value x (i) ; 3. for each component, approximate the relation between its sampled coefficient and the corresponding sampled model inputs by means of a GP surrogate model;

4. formulate the POD+GP surrogate model ĥ(x) as the weighted sum of the more significant POD components where weights are the GP surrogate models depending on x.

Reduction of the output dimension by proper orthogonal decomposition (POD)

The key idea of the snapshot method [START_REF] Siade | Water resources research. Snapshot selection for groundwater model reduction using proper orthogonal decomposition[END_REF] is to achieve a POD of the centred snapshot matrix

H = h (i) j -N -1 l N l k=1 h (k) j 1≤i≤N l 1≤j≤p ∈ M N l ,p (
R), which gathers the water level computed at each mesh point for the N l snapshots, from which the sample mean is substracted.

Based on many observations of a random vector, the POD gives the orthogonal directions of largest variances (or modes) in the probabilistic vector space in order to reduce the vector space dimension [START_REF] Chatterjee | Current Science. An introduction to the proper orthogonal decomposition[END_REF]. Note that for simplicity purpose, the adjective centred is dropped in the following when referring to the centred snapshot matrix H.

The POD of the snapshot covariance matrix

C = N -1 l H T H ∈ M p (R) is equivalent to the Singular Value Decomposition (SVD) of the snapshot matrix H: H = U Λ V T = rp k=1 λ k u k v T k , (4)
where U ∈ M N l (R) is an orthogonal matrix diagonalizing HH T (u k , the k th column of U, is a left singular vector of H), where V ∈ M p (R) is an orthogonal matrix diagonalizing H T H (v k , the k th column of V, is a right singular vector of H), and where Λ ∈ M N l ,p (R) is a rectangular diagonal matrix including r p = min(N l , p) singular values on its diagonal. The singular values {λ k } 1≤k≤rp are the square roots of the eigenvalues of C. Note that in this study, since the size of the training set N l is lower than the number of mesh points p = 41 416, the rank of H is here r p = N l . At the k th mesh point, the snapshot h k (x (i) ) can then be retrieved as a linear combination of r p modes {Ψ i } 1≤i≤rp :

h k x (i) = U Λ V T ki = U k: Λ V T :i = rp j=1 γ k,j Ψ j x (i) , (5)
where for any j ∈ {1, . . . , N l }, γ p,j := U k,j and Ψ j x (i) := ΛV T j,i . From that, we want to approximate each relation x → Ψ j (x) by a GP surrogate model Ψ gp,j from the dataset x (i) , Ψ j x (i) 1≤i≤N l in order to propose the following POD+GP surrogate model:

h k (x) = rp i=1 γ k,i Ψ gp,i (x), (6)
This POD+GP surrogate model requires the construction of r p GP surrogate models.

Learning of the significant POD modes by Gaussian process (GP) modelling

As stated by [START_REF] Rasmussen | Gaussian Processes for Machine Learning[END_REF], a GP is a random process (here the mode Ψ i ) indexed over a domain (here R d ), for which any finite collection of process values (here Ψ i x (j) 1≤j≤N l ) has a joint Gaussian distribution. Concretely, let Ψ i be a Gaussian random process fully described by its zero mean and its correlation π i :

Ψ i (x) ∼ GP 0, σ 2 i π i (x, x ) , (7)
with π i (x, x ) = E Ψ i (x) Ψ i (x ) . In our case, the correlation function π (or kernel) is chosen as a squared exponential:

π i (x, x ) = exp - x -x 2 2 2 i , ( 8
)
where i is a length scale describing dependencies of model output between two input vectors x and x , and where σ 2 i is the variance of the output signal. Squared exponential kernel leads to satisfying results but other kernel functions could have been considered, such as a decreasing exponential one or a Matérn one -with their associated hyper-parameters. The choice of the kernel is still an open problem and can be mitigated using the available information on the problem. The square exponential kernel leads to very smooth, thus stable results. Furthermore, it implies that the model is exact at sample points; it does not introduce any other strong assumptions, hence its wide usage among practitioners.

Then the surrogate model of interest is the mean of the GP resulting of conditioning Ψ i by the training set

Ψ i x (k) 1≤k≤N l . For any x * ∈ R d , Ψ gp,i (x) = N k=1 β k,i π i x, x (k) , (9)
where

β k,i = (Π i + τ 2 I N l ) -1 Ψ i x (1) . . . Ψ i x (N l ) T with Π i = π i x (j) , x (k) 1≤j,k≤N l
, and where τ (referred to as the nugget effect) avoids ill-conditioning issues for the matrix Π. The hyperparameters { i , σ i , τ } are optimized by maximum likelihood applied to the data set D N using the L-BFGS-B algorithm [START_REF] Wales | The Journal of Physical Chemistry A. Global Optimization by Basin-Hopping and the Lowest Energy Structures of Lennard-Jones Clusters Containing up to 110[END_REF].

Quality measures for the POD+GP surrogate model

In the present study, two common error metrics are used to assess the quality of the surrogate water level both on the entire mesh (global approach) and at each point of the mesh (local approach): the root mean square error (RMSE) and the predictive coefficient (Q 2 ). This validation is carried out over an input-output validation dataset D v of size N v .

Root mean square error (RMSE)

The RMSE is used to measure the accuracy of the model, it should be 0 when the model is perfect. At the k th given mesh node, it is defined as the square root of the mean square error (MSE) measuring the square distance between the surrogate model and the reference model:

MSE k (D v ) = N -1 v Nv i=1 h (i) k -h (i) k 2 and RMSE k (D v ) = MSE k (D v ) (10)
Their global counterpart are:

MSE(D v ) = p -1 p k=1 MSE k (D v ) and RMSE(D v ) = MSE(D v ). Predictive coefficient (Q 2 )
At the k th mesh node, the Q 2 predictive coefficient is defined as:

Q 2,k = 1 - MSE k (D v ) MSE k (D v ; mean) ( 11
)
where

MSE k (D v ; mean) = N -1 v Nv i=1 h (i) k -h (i) 2
is the MSE of the "averaging model" returning the mean of the learning outputs whatever the input parameter value.

The global counterpart of MSE

k (D v ; mean) is MSE(D v ; mean) = p -1 p k=1 MSE k (D v ; mean). Thus, the global counterpart of Q 2,k is: Q 2 = 1 - MSE(D v ) MSE(D v ; mean) . ( 12
)
The predictive coefficient measures the performance of the surrogate model with respect to the simplest one which consists in averaging the learning output values. When Q 2 is lower than (resp. equal to) zero, the surrogate is worse than (resp. equal to) the learning output values average. When Q 2 is equal to one, the surrogate interpolates the validation dataset. In practice, the surrogate is deemed appropriate when Q 2 is greater than 0.8. The predictive coefficient is also found under the name of Nash-Sutcliffe model efficiency coefficient in the hydrological literature where is assesses the predictive capacity of the simulated discharge over a time window with respect to observed discharges [START_REF] Nash | River flow forecasting through conceptual models part I: A discussion of principles[END_REF].

Quantify and explain the output uncertainty due to input uncertainty propagation

Once the model is built and validated, it can be used instead of the reference model in an uncertainty quantification study. After propagating the input uncertainties through the surrogate model by means of specific Monte Carlo methods, we can conduct a statistical analysis on the output uncertainty as well as a sensitivity analysis to explain how the uncertain input parameters contribute to this output variability.

Statistical analysis on the output

Using a standard MC approach on the validation data set D v , the mean value and standard deviation of the water level at the k th mesh point are formulated as:

µ k = 1 N v Nv i=1 ĥ(i) k and σ k = 1 N v -1 Nv i=1 ĥ(i) k -µ k 2 . ( 13
)

Sensitivity analysis on the output with respect to the inputs

Sobol' indices [START_REF] Sobol | [END_REF] are commonly used for sensitivity analysis. They provide the shares of the QoI variance V attributable to the different uncertain inputs. Under the hypotheses that random input variables are independent, here the roughness coefficients and the upstream flow, and the random QoI is square integrable, here the water level h, the decomposition of the QoI reads:

V = d i=1 V {i} + d j=i+1 V {i,j} + • • • + V {1,2,...,d} = J⊂{1,2,...,d} V J , (14)
where

V := Var [QoI], V i := V [E[QoI|X i )], V ij := V [E[QoI|X i X j ]] -V i -V j
and more generally, for any I ⊂ {1, . . . , d},

V I := V [E[QoI|x I ]] -J⊂I s.t. J =I V J .
Then, we obtain:

1 = d i=1 S {i} + d j=i+1 S {i,j} + • • • + S {1,2,...,d} = J⊂{1,2,...,d} S J , (15)
where for any J ⊂ {1, 2, . . . , d}, S J = V J V is called a Sobol' index and belongs to the interval [0, 1]. S {i } is the first order Sobol' index corresponding to the ratio of output variance due to the i th input parameter uniquely, and S {ij} is the second-order Sobol' index describing the ratio of output variance due to the i th parameter in interaction with the j th parameter. Also the total Sobol' index that corresponds to the whole contribution of the i th input parameter reads:

S T i = I⊂{1,...,d} I i S I . (16)
The computation of first order Sobol' indices requires simple integration, those of the second order requires double integration, and so on. Many Monte Carlo techniques exist to estimate these integrals. In this study, the Sobol' indices are estimated using the algorithm proposed in [START_REF] Saltelli | Computer Physics Communications. Variance based sensitivity analysis of model output. Design and estimator for the total sensitivity index[END_REF].

Lastly, we note that the expression ( 14) is defined for a scalar QoI, such has the water level h k at mesh node k ∈ {1, 2, . . . , p} where p is the mesh size. Consequently, we can easily plot the different Sobol' indices over the mesh on which is defined the model output. Furthermore, this information can be summarized using the generalized Sobol' indices [START_REF] Lamboni | Reliability Engineering & System Safety. Multivariate sensitivity analysis to measure global contribution of input factors in dynamic models[END_REF]:

∀J ⊂ {1, 2, . . . , d}, S J = p k=1 V [k] S [k],J p =1 V ( ) [k] . ( 17
)
4 RESULTS The design of experiment (DoE) for the training and validation data set was generated using Latin Hypercube Sampling (LHS) [START_REF] Mckay | Technometrics. A comparison of three methods for selecting values of input variables in the analysis of output from a computer code[END_REF] which is a statistical method for generating a nearrandom sample of parameter values from a multidimensional distribution. Considering d the number of input variables, LHS strategy scales as o(d) while other strategies require a larger number of samples; for instance to insure the convergence of first order statistics [START_REF] Damblin | Journal of simulation. Numerical studies of space filling designs: optimization of Latin Hypercube Samples and subprojection properties[END_REF]. The LHS space-filling experimental design is shown in Fig. 2, it is associated with a limited computational cost. While more advanced sampling method could be used, LHS strategy was deemed efficient for the present study.

Learning and test samples

Surrogate model

The LHS strategy has been applied twice. A first time to build a 2000-sample training set and a second one to create a 1000-sample validation set. Here, the GP kernel was prescribed to a Matern(2.5) function. The validation set was only used to assess the quality of the surrogate model with RMSE and Q 2 error metrics. The dimension of the output space was reduced with a POD in order to limit the cost of the GP surrogate. The number of modes, also called principal components, to be retained is justified by two criteria taken into account:

• Elbow criterion: on the scree of the POD modes, there is a decrease (elbow) followed by a more regular decrease. In our case, as shown in figure 3, a decrease occurs at the fourth mode, then a regular decrease from the fifth mode. Thus only the first four modes are of interest.

• Kaiser's criterion: only those modes whose inertia is greater than the average inertia should be retained. This criterion leads us to select 4 modes, explaining 99.6% of the total inertia. Indeed, the first principal component explains 95.84% of the total inertia, the second 1.80%, the third 1.46% and the fourth 0.48%.

studyThe cost of the GP surrogate significantly decreases when the dimension of the output is reduced5 (6 times smaller) applying the POD, as presented in Tab. 2. The output dimension is indeed reduced from 41 416 elements to 4 components that explain 99.6% of the variance of the QoI. But, the physical interpretation of the different modes is not always obvious. As displayed in Fig. 4, the first mode, which explains 95.84% of the output variability, seems to represent the effect of the upstream discharge on the average water level height. Indeed, this component is essentially negative thus its weighting will increase everywhere the average water level height if negative coefficient or decrease everywhere the average water level if positive coefficient. While the second mode, which explains 1.80% of the output variability, seems to The POD+GP surrogate quality is very good with respect to global error metrics RMSE = 0.8 cm and Q 2 = 0.99748. Locally, the quality deteriorates near the boundary of the catchment area as well as along the dikes as shown in Fig. 5. The heterogeneity of the mesh with small cells in the river bed (≤ 40 m), near the dikes (≤ 80 m) and larger cells in the flood plain (≤ 150 m) should be noted and may hide some local failures of the surrogate in the global RMSE and Q 2 criteria.

Sensitivity analysis with the POD+GP surrogate

The POD+GP surrogate is used to carry out a variance-based sensitivity analysis (SA) over the entire simulated area, with a focus on mesh node 29 515 where Marmande, a city prone to flooding, is located. The POD+GP surrogate allows for a reliable estimation of first and second order statistical moments at Marmande as shown in Tab. 3: the water level mean and standard deviation estimated from the direct model T2D and surrogate are in good agreement with an under estimation of 1.3% for the mean computed with the surrogate.

Given the statistical distributions for the input variables, the SA at Marmande highlights that most of the water level variance is explained by the upstream discharge Q and to a lesser extend, by the Strickler friction coefficient K s,4 prescribed between Marmande and La Role as displayed in Fig. 6. At this location, the floodplain friction coefficient K s,1 and the friction coefficients upstream of Mas d'Agenais (K s,2 ) and upstream of Marmande (K s,3 ) have barely no impact on the water level. It should be noted that the bootstrap method [START_REF] Archer | Sensitivity measures, anova-like techniques and the use of bootstrap[END_REF] is used to estimate the variance of the Sobol' indices, this variance is represented by the black error bars in Fig. 6. These indicate that the computation of the SA indices is converged and reliable. It should also be noted that for each input variable, the first (S) and total Sobol' (S T ) indices at Marmande are equal, meaning that, at this location, the multivariate impact of the input on the water level is minimal. The 2D computation and representation of first Sobol' indices confirms that the variance of the water level is mostly explained (81 % on average) by the input discharge as shown in Fig. 8 and Fig. 9. The floodplain friction coefficient K s,1 has no impact on the analysis, the upstream friction coefficient in the river bed K s,2 has a small impact on the water level close to Mas d'Agenais, the friction coefficient K s,3 between Mas d'Agenais and Marmande explains up to 10 % of the water level variance close to Marmande and the downstream friction coefficient K s,4 has an impact over the entire domain with most significance at the upstream and downstream boundaries. As the sum of the first Sobol' indices is smaller than 1, higher order Sobol' indices are non zero, meaning that multivariate effects between Q and K s explain the remaining part of the water level variance.

CONCLUSION AND PERSPECTIVES

In this paper, an uncertainty quantification study was carried out with a 2D numerical solver for the Shallow Water Equations on a section of the Garonne river. It consisted in building a Gaussian process surrogate model on a POD-reduced 2D water level output field.

The surrogate model was formulated with respect to friction coefficients and input discharge, the distribution for friction is supposed to be uniform and centered around calibration values while the discharge distribution is supposed to be Gaussian, centered around a high flood value. The construction of the surrogate was achieved over a 2000-sample training data set and it was validated over a 1000-sample data set. The dimension of the quantity of interest was reduced from 41416 elements to 4 principal components using the POD which has resulted in a significant reduction of the computational cost of the surrogate. The correlation kernel was here prescribed as a Matern(2.5) function. The quality of the POD+GP surrogate model was assessed, the surrogate was deemed satisfying with Q 2 metrics close to 1 for the entire domain and RMSE smaller than 0.01m. The quality of the surrogate decreases near the dikes. The surrogate was used to perform a global sensitivity analysis based on variance decomposition. It was demonstrated that the upstream discharge is the predominant input variable and explains more than 80 % of the water level variance. The downstream friction coefficient is also a significant input with heterogeneous influence.

It is essential to mention that the conclusions for this study are strongly related to the hypothesis made for the statistical distribution of the inputs. For instance, further study should investigate wider ranges for flood plain coefficients that are highly unknown and may significantly over time as flood events occur on the catchment.

The results of the sensitivity analysis allows for a better understanding of the physics as well as classification of major sources of uncertainty. The latter is of great importance in the context of data assimilation where the control vector should be properly defined to include key factor to improve the model outputs. It was here highlighted that in order to improve water level at Marmande, the control vector should include at least the upstream discharge and the downstream friction coefficient. A perspective for this study thus stands in the implementation of an ensemble-based data assimilation algorithm to improve input discharge and friction assimilation water level observations in the system. Additionally, the cost of the ensemble integration should be reduced using the surrogate model in place of the direct hydraulic solver.

1 3

 1 w /ρ air [kg.m -3 ] is the water/air density, P atm [Pa] is the atmospheric pressure, U w,x and U w,y [m.s -1 ] are the horizontal wind velocity components, C DZ [-] is the wind influence coefficient, K s [m .s -1 ] is the river bed and floodplain friction coefficient, using the Strickler formulation [?]. F x and F y [m.s -2 ] are the horizontal components of external forces (friction, wind and atmospheric forces), h [m] is the water depth, H [m] is the water level

Figure 1 :

 1 Figure 1: Study area of the Garonne river

Figure 2 :

 2 Figure 2: Latin Hypercube Sampling (LHS) DoE for a 300-sample data set, along (K s,1 , K s,3 , Q) directions on the left panel and along (K s,3 , Q) directions on the right panel.

Figure 3 :

 3 Figure 3: Principal component analysis

Figure 4 :

 4 Figure 4: Principal component analysis

Figure 5 :

 5 Figure 5: Root mean square error

Fig. 7

 7 Fig.7displays the mean and the standard deviation of the water level over the 2D domain estimated with the POD+GP surrogate. The mean varies between 0 m near the limits of the domain and 21.57 m at Marmande, where it reaches its maximum. In the floodplain, the mean water level ranges from 3.21 m and 7.8m close to the dikes. The water level standard deviation ranges from 0 m to 0.4 m in the floodplain between Mas d'Agenais and Marmande, where the flow is highly bi-dimensional.

Figure 6 :

 6 Figure 6: First (S) and total (S T ) Sobol' indices estimated with POD+GP surrogate at Marmande.

Figure 7 :

 7 Figure 7: a-Mean water level, b-standard deviation of the water level, estimated with the POD+GP surrogate.

Figure 8 :

 8 Figure 8: First order Sobol' sensitivity indices computed with the POD+GP surrogate with respect to the input discharge Q.

Figure 9 :

 9 Figure 9: First order Sobol' sensitivity indices computed with the POD-GP surrogate with respect to the friction coefficients K s,1 (flood plain), K s,2 upstream Mas d'Agenais, K s,3 (between Mas d'Agenais and Marmande and K s,4 (downstream of Marmande).

Table 2 :

 2 CPU run time comparison between GP without and with POD.

Table 3 :

 3 Statistical moments of the water level height in Marmande.

More information can be found on the website www.opentelemac.org.

BATMAN-OT can be downloaded from https://gitlab.com/cerfacs/batman.

More information on: http://www.openturns.org.

More information on: https://scikit-learn.org.

The remaining cost can be considered significant compared to linear surrogate models such as polynomial chaos expansion. This situation is well-known and naturally explained by the learning sample size increasing the cost of inverting the covariance matrix. For prediction, this surrogate model is as fast as the others and also provides a measure of its error.

ACKNOWLEDGEMENTS

We would like to thank Nicole Goutal from EDF R&D's LNHE department for providing us with the case study, Pamphile Roy from CERFACS for his algorithmic developments in the BATMAN-OpenTURNS library and the EDF R&D's PRISME department for their support in using OpenTURNS. We thank all these people for the enriching discussions we were able to have on uncertainty methods for two-dimensional river hydraulics issues.


Rachid Ouaret
email: rachid.ouaret@toulouse-inp.fr

Pascal Floquet

Jean-Pierre Belaud

Stéphane Negny

Interval-Valued Data Learning for Robustness of Energy Recovery Systems

Keywords: Machine Learning for Symbolic Data, Interval-valued data, Flexibility and Robustness, Heat Exchanger Networks (HEN)

Modeling of energy recovery systems, Heat Exchanger Networks (HEN) as an example, are complex processes due to many parameters involved and the lack of knowledge on the impact of operational variables on the response. These variables are often affected by different types of uncertainty as to measurement errors, computation errors, or imprecision related to the underlying method. The uncertainty in the data may be treated by considering, rather than a single value for each variable, the interval of values in which it may fall, histograms, or other multivalued data: symbolic data. This work aims to use, when possible, the symbolic data analysis to adapt the classical mathematical HEN models. It deals with the study of continuous interval data through suitable Principal Component Analyses and Regression for two purposes: clustering exchanges (i) classification of exchangers to detect those impacted by uncertainty factors and (ii) evaluation of the relationship between the different process parameters (inlet temperature, heat transfer coefficient, etc.) on interval data. The new method has been tested on a real data set and the numerical results are reported. The symbolic approach provides a simple way to study a great number of scenarios.

Introduction

Energy saving is an important issue for both industries and society. In the industrial chemical process, Heat Exchanger Networks (HEN) are widely used techniques for reducing external heating and cooling utilities. Data generated by those complex systems has increased drastically over the past few years. Suppose we have 20 exchangers on which 4 variables are measured (2 input and 2 output hot and cold input streams). If we assume that we have hourly values of these variables for 1 year, then each exchanger is described by 35040 data. If, on the other hand, certain characteristics of the exchangers (let's say 3) vary over time, then a tabular representation conventionally used in data analysis would contain 2 102 400 values (35040 × 20 × 3). This is not very large compared to what could be provided by the industry in real-time.

Even such data is ubiquitous for a larger scale of HEN systems, data-oriented based approaches are often analyzed with simplified models by aggregating data. During this operation, some information on the variability aspects is lost. The robustness assessment of HEN is therefore affected since its flexibility is conditioned by the variability of the uncertain parameters. When the size n of entries (exchangers) and p number of variables (features) are very large, classical analysis can be problematic.

To address this problem, we use a more complex object description to capture the variability of measured parameters on each exchanger. When dealing with quantitative variables, complete information can be achieved by describing a set of statistical units in terms of interval data, histogram, ... rather than a single-valued variable. Mathematically, intervalvalued data with measurements on p random variables, are p-dimensional hyperrectangles in R p . Such data need to be visualized, synthesized, and compared on factor spaces. This paper is in a complementary perspective to Floquet et al's work [START_REF] Floquet | Analysis of operational heat exchanger network robustness via interval arithmetic[END_REF]. They have initiated the robustness analysis of a simple exchange networks using interval arithmetic. They pointed out the "butterfly effect" of the alteration of characteristics of some heat exchangers on the operation of the HEN. Indeed, a maintenance operation of an exchanger, alter the variability of the other parameters of the network. Reducing fouling would imply variations in pressure and shift flows between parallel branches, all changing over time in a way that is difficult to predict [START_REF] Macchietto | Energy recovery in heat exchanger networks in a dynamic, bigdata world: Design, monitoring, diagnosis and operation[END_REF]. The purpose of this paper is to check whether there are groups of exchangers characterized by the same properties in terms of their responses to external fluctuations. Following the same idea (interval-valued data), the Symbolic Data Analysis (SDA) [START_REF] Billard | Symbolic Data Analysis: Conceptual Statistics and Data Mining[END_REF][START_REF] Bock | Analysis of symbolic data: exploratory methods for extracting statistical information from complex data[END_REF] is used to study how uncertainty in the output of a model can be apportioned to different sources of uncertainty.

Machine Learning for symbolic data

Figure 1 outlines the design and methodological scheme of the proposed method. This methodology draws on two main components: (B) Symbolic data through interval arithmetic, and (C) machine learning for symbolic data. Part A corresponds to the HEN model to be studied or simulated. The interval-valued data is then constructed in part (B).The core function of the proposed method included in this paper is part (C). The originality of this research lies in the combination of a traditional robustness analysis of HEN with SDA.

Figure 1: Overview of the SDA integration and data-driven processing for HEN robustness analysis.

Symbolic Data and statistics of interval-valued variables

In classical statistics, n × p data matrix X = X i j is defined between n individuals and p variables, where each cell (i, j) contains a unique value x i j . A symbolic objects are more complex than a simple valued variable description, symbolic data can contain internal variation of the features representing imprecise knoweledge and can be structured. The symbolic analysis generalizes the classical data analysis, e.g. To illustrate this point with a simple HEN example, if the outlet temperature of two exchangers ω 1 , ω 2 is given by T out (ω 1 ) = 90

• C, T out (ω 2 ) = 120 Simple statistical (mean, variance, ...) for interval-valued variables have been defined in [START_REF] Bertrand | Descriptive statistics for symbolic data[END_REF]. Let consider Y j ≡ Z be the j th intervalvalued random variable, and

Z (ω u ) = [a u , b u ] is a realization of Z for the observation ω u over the observed interval [a u , b u ]. The empirical distribution function, F Z (ξ ), is the distribution function of a mixture of m distributions {Z (ω u ) , u = 1, 2, . . . , m}.
The central and dispersion parameters of a variable all derived from a strong assumption: the inherent fluctuation within random intervals and rectangles is uniformly distributed:

f (ξ ) = 1 m ∑ u:ξ ∈Z(u) 1 (b u -a u )
. The symbolic sample mean for interval-valued data is given by

Z = 1 2m ∑ u (b u + a u ) , (1)
and the sample variance is given by

S 2 = 1 3m ∑ u b 2 u + ba u + a 2 u - 1 4m 2 ∑ u (b u + a u ) 2 .
(2)

PCA for interval-valued data

A principal component analysis is designed to reduce p-dimensional observations into s-dimensional components (where s p) in an interpretable way, such that most of the information in the data is preserved. Let i = 1, . . . , n denote n objects (exchangers) described by p features (or variable) Y 1 , . . . ,Y p (temperatures, ...). The symbolic data matrix used for interval PCA is given by X

=    ξ 11 • • • ξ 11 . . . . . . . . . ξ n1 • • • ξ np   
, where ξ i j = x i j , x i j is the interval of possible values of variable j for the exchager i, and the symbolic data vector can be denoted by x x x i = (ξ i1 , . . . , ξ ip ) = x i1 , x i1 , . . . , x ip , x ip . The data point is represented in R p space by hyperrectangles R i with 2 p vertices. There are mainly two methods to solve the algebraic mapping to lower dimension: Vertices and Centers methods. To find the factorial axes for Centers method, a classical PCA is applied to the centers c i ∈ R p of the n hyperrectangles R i . The coordinates of i th center c i is denoted by x c i j , where x c i j is computed for i = 1, . . . , n and j = 1, . . . , p. From that case, the centers of hyperrectangles in its rows is an n × p data matrix, denoted by X and the j th column of X is denoted by the feature Y c j . The interval principal components values are obtained by computing first the υ th principal component of the center c i is given by

ψ c iυ = p ∑ j=1 x c i j -x c i j • η jυ (3)
where η jυ = (η 1υ , . . . , η pυ ) is the υ th eigenvector of S (sample covariance matrix associated with the dataset). It is possible to find an ψ iυ , ψ iυ of the possible values of υ th principal component ψ c iυ of c i [START_REF] Cazes | Extension de l'analyse en composantes principales à des données de type intervalle[END_REF]. For the object i :

ψ iυ = p ∑ j=1 max x i j ≤x r i j ≤x i j
x r i jx c j • η jυ (4)

ψ iυ = p ∑ j=1 min x i j ≤x r i j ≤x i j x r i j -x c j • η jυ (5)
3. Application to HEN data

Data simulation design and preliminary analysis

The starting point for this application is the data table from step B in Figure 1, which is the result of HEN model analysis by interval arithmetic. In order to evaluate the behavior of the exchanges on a more complex network, simulated data on the initial table were undertaken and the simulation procedure is as follows:

1. For each variable (T h in min h in , max h in , T c in [min c in , max c in ], T h out min h out , max h out , T c out [min c out , max c out ]
) and for all exchagers, we look for the minimum and the maximum.

2. Simulate 40 points (exchangers) using random variable from the uniform distribution. For example, for the exchanger No. 30 (E30), the minimum input temperature for the cold stream is obtained by the following probabilistic simulation scheme: T c in (E30) ∼ U (min (min c in ) , max (min c in )). Retrieve the interval matrix of the 4 variables (T h in ,T c in ,T h out ,T c out ).

Statistical description for the interval-valued data of the initial matrix and the simulated one, using equations 1 and 2, is presented in the following Table 1:

Table 1: Descriptive statistics (mean and standard deviation s.d) for the interval-valued data of HEN Mean [124.20:131.42] [46.32:61.30] [70.32:84.30] [84.68:94.95] sd [44.47:45.42] [38.25:23.67] [20.81:12.03] [47.34:38.20] Simulation Mean [112.14:134.02] [45.34:42.93] [68.95:81.42] [76.98:97.71] sd [20.77:28.33] [22.91:23.85] [11.38:9.02] [33.34:25.32]

T h in [T c in ] T h out [T c out ] Initial data

Principal Component Analysis for interval-valued data

The interpretation of the position of the interval-valued data in the principal plane is the same as in the classical principal component analysis situation. In Figure 2, we show the results with respect to the first three axes, achieved by the Symbolic PCA using 4 and 5 (centers method). Notice that the 61% of the total inertia is explained by the first two axes in the case of simulation (40 exchangers), and 98.5% of in the case of initial data (4 exchangers). In Figure 2, closeness among clusters exchanger mainly influenced by the same descriptors. Only one group consisting of exchangers 1 and 3 can be identified. The constituent elements of this cluster are influenced by the same main factors. Exchangers 2 and 4 are detached from the cluster from a 3D perspective. This observation is in line with the conclusions of the initial study [START_REF] Floquet | Analysis of operational heat exchanger network robustness via interval arithmetic[END_REF].

For simulated data, the same analysis can be achieved using additional information, and at this stage, it is also difficult to give an interpretation of the similarity in size and shape among exchangers. In order to shed new light on the variability of outlet temperatures as a function of input streams, we propose a linear regression on interval values.

Regression and prediction for interval-valued data

The model proposed here is based on the classical formation of HEN models. However, a notable difference lies in the integration of symbolic objects in the model. To understand how the input (cold and hot) temperature stream affects the output temperature, we propose a regression model for interval-valued temperature. For a hot stream:

T h out min h out , max h out = β 0 + β 1 • T h in min h in , max h in + β 2 • T c in [min c out , max c out ] + ε (6)
The following results are obtained using the simulated data in the Table 2. which is a set of the rest 5 exchangers used to provide an unbiased evaluation of the estimated model. The Table 3 shows the results of the predictions using the fitted model and corresponding data test. The size of the coefficient for each independent variable (T h in min h in , max h in and T c out [min c out , max c out ]) gives the size of the effect that variable is having on the dependent variable T h out and T c out . The estimated coefficients tells how much the output stream is expected to increase when that inputs streams (hot and cold) increases by one. For example the coefficients for the model 6, T h in differed by 1 • C (and T c in did not differ) T c out will differ by 0.17 • C units, on average. The estimated intercept β 0 , is the expected mean value of T h out and T c out when all inputs are 0. In our experiments, the T h in and T c out never comes close to 0, then intercept has no meaningful interpretation.

Conclusion

Symbolic Data Analysis (SDA) extends statistics and multivariate data analysis to deal with data structured in a distributional form with complex internal variations. In this paper, comprehensive modeling via SDA has been presented that moves substantially beyond the traditional modeling in HEN robustness analysis. It proposed some new approaches which intended to redefine the robustness of HEN based on interval data. First, the detection of Exchangers cluster which would be affected by common factors has been modeled by Symbolic Principal Component Analysis. Second, the relation between inherent variation, expressed by intervals, of input and output temperature has been modeled using the linear regression method for interval-valued variables. Future research should include histogram valued-data.

  x = c, c ∈ R is equivalent to the symbolic interval ξ = [c, c]. A full conceptualization of symbolic objects can be found in Bock and Diday (2012). Let Ω be a set of individuals, D containing the descriptions of individuals and the descriptions of classes of individuals, a, a mapping defined from Ω into D which associates to each ω ∈ Ω a description d ∈ D by using intervals, histograms, etc. More formally, symbolic object is a triplet s = (a, R, d) where R is a relation between descriptions, d is a description, and a is a mapping defined from Ω in L depending on R Interval-Valued Data Learning for Robustness of Energy Recovery Systems 3 and d. For instance L = {true, f alse} or L = [0, 1] and R may be one of the relations in {=, ≡, ≤, ⊆, =⇒}.

  • C the description of the class D = (ω 1 , ω 2 ) obtained by a generalization process can be [90, 120]. The symbolic object s is defined by a triple s = (a, R, d) where d = [90, 120], R = ∈ and a is the mapping: Ω → {true, f alse} such that a is the true value of T out (ω) Rd,written a (ω) = T out (ω) ∈ [90, 120]. An individual ω is in the extent of s if and only if a (ω) = true.

Figure 2 :

 2 Figure 2: Principal 3D-space with data of interval type of HEN. Factorial for 4 Exchangers (left) and for 40 Exchangers (right).

Table 2

 2

		: Results of interval-valued data regression
		Dependent variable
		T h out min h out , max h out	T c out [min c out , max c out ]
	β 0	63.22	37.46
	β 1	0.0046	0.17
	β 2	0.244	0.6
	To evaluate this approach, two datasets were used: (i) training dataset which is a set
	of the first 35 exchangers used to fit the parameters of the model 6 and (ii) test dataset

Table 3 :

 3 The prediction results Predictions Test .01] [61.41:80.81] [73.61:67.64] [39.89:119.97] E37 [73.25:82.21] [76.12:104.68] [69.72:89.61] [115.31:89.80] E38 [68.06:76.77] [63.75:87.63] [57.47:91.65] [46.96:104.48] E39 [64.73:77.73] [69.71:80.77] [88.25:85.72] [113.54:66.64] E40 [60.93:72.54] [52.76:75.85] [52.70:88.04] [83.85:100.02]

	Exchangers	T h out	T c out	T h out	[T c out ]
	E36	[64.21:75