Cheng-Wei Chou

Ping-Chiang Chou

Jean-Joseph Christophe

Adrien Couetoux

Pierre De Freminville

Nicolas Galichet

Chang-Shing Lee

Jialin Liu

David Saint-Pierre

Michèle Sebag

Strategic Choices in Optimization

à 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.

Strategic Choices in Optimization

1. INTRODUCTION


A Meilland

Ph Stee

O Chesneau

C Jones

VLTI/MIDI observations of 7 classical Be stars ⋆

Keywords: high angular resolution -Techniques: interferometric -Stars: emission-line, Be -Stars: winds, outflows, disks -Stars: individual (p Car, ζ Tau, κ CMa, α Col, δ Cen, β CMi, α Ara

Context. Classical Be stars are hot non-supergiant stars surrounded by a gaseous circumstellar envelope that is responsible for the observed IR-excess and emission lines. The origin of this envelope, its geometry, and kinematics have been debated for some time. Aims. We measured the mid-infrared extension of the gaseous disk surrounding seven of the closest Be stars in order to constrain the geometry of their circumstellar environments and to try to infer physical parameters characterizing these disks. Methods. Long baseline interferometry is the only technique that enables spatial resolution of the circumstellar environment of classical Be stars. We used the VLTI/MIDI instrument with baselines up to 130 m to obtain an angular resolution of about 15 mas in the N band and compared our results with previous K band measurements obtained with the VLTI/AMBER instrument and/or the CHARA interferometer.

Results. We obtained one calibrated visibility measurement for each of the four stars, p Car, ζ Tau, κ CMa, and α Col, two for δ Cen and β CMi, and three for α Ara. Almost all targets remain unresolved even with the largest VLTI baseline of 130m, evidence that their circumstellar disk extension is less than 10 mas. The only exception is α Ara, which is clearly resolved and well-fitted by an elliptical envelope with a major axis a=5.8±0.8 mas and an axis ratio a/b=2.4±1 at 8 µm. This extension is similar to the size and flattening measured with the VLTI/AMBER instrument in the K band at 2 µm. Conclusions. The size of the circumstellar envelopes for these classical Be stars does not seem to vary strongly on the observed wavelength between 8 and 12 µm. Moreover, the size and shape of α Ara's disk is almost identical at 2, 8, and 12 µm. We expected that longer wavelengths probe cooler regions and correspondingly larger envelopes, but this is clearly not the case from these measurements. For α Ara this could come from to disk truncation by a small companion ; however, other explanations are needed for the other targets.

Introduction

Classical Be stars are hot non-supergiant stars that are surrounded by a gaseous circumstellar environment responsible for the presence of an IR-excess and emission lines in their spectra. A generally accepted view of this material is a dense equatorial region dominated by rotation, often called a "disk" where most of the IR-excess and emission lines originate, surrounded by a more diluted, radiatively driven, polar wind. Despite the fact that these systems have been studied for decades, the disk and wind formation, geometry, and kinematics, remain contentious issues.

First VLTI observations of 4 Be stars, namely Achernar (Domiciano de Souza et al. 2003), α Ara (Meilland et al. 2007a), κ CMa (Meilland et al. 2007b), and δ Cen (Meilland et al. 2008) have shown evidence of different geometries likely to due different mass-ejection processes such as fast rotation, radiatively driven winds or binarity. Thus, it is not clear if Be stars can be considered as a homogeneous group of stars in terms of massejection processes (Stee & Meilland 2009).

Measuring the N-band extension of the circumstellar environments of classical Be stars can put strong constraints on their physical parameters as shown by Chesneau et al. (2005) for α Ara. For this star, VLTI/MIDI observations with 100 m base-Send offprint requests to: meilland@mpifr-bonn.mpg.de ⋆ Based on observations made with ESO Telescopes at Paranal under programs 073.A-9014G, 077.D-0097 and 082.D-0189 lines led to nearly unresolved visibilities, demonstrating that the envelope was even smaller than what was previously predicted by Stee & Bittar (2001). These authors used long term spectroscopic variability measurements and concluded that the disk is likely truncated by a companion.

We initiated an observational campaign on the brightest, closest, southern hemisphere classical Be stars using the VLTI instruments MIDI and AMBER. In this paper we present new spectrally resolved mid-infrared VLTI/MIDI interferometric measurements for seven classical Be stars, p Car, ζ Tau, β CMi, κ CMa, α Col, δ Cen, and α Ara.

The paper is organized as follows. In Sect. 2 we present the observations and the data reduction process. In Sect. 3, the calibrated visibilities are presented and then compared with a model composed of a unresolved star plus uniform disk. Our results and conclusions are discussed in Sect. 4, and Sect. 5, respectively. HD107446, HD151249, HD163376, HD169916, HD198048, HD217902, and HD218594. The observation log is presented in Table 1.

Observations and data reduction

The 2004 observations of δ Cen were obtained with UT1-UT3 in HIGH-SENS mode. All of the other observations were taken with the UT1-UT4 baseline using the SCI-PHOT mode to enable a better visibility calibration (Chesneau 2007). The PRISM low-spectral dispersion mode, allowed us to obtain spectrally resolved visibilities with R=30 across the N band. We have used the standard MIA + EWS, (version 1.6) package to reduce the data and estimate the errors bars (see Fig. 1).

Results

The calibrated visibilities

The N band calibrated visibilities for the program Be stars are plotted as a function of wavelength in Fig 1 . The only clearly resolved (i.e., V+3σ <1) target is α Ara along the B 1 baseline. It is also barely resolved (i.e., V+σ <1) along the other two baselines. Moreover, δ Cen (for both baselines) and ζ Tau (for λ >10µm) are also barely resolved. All of the other targets, κ CMa, p Car, α Col, and β CMi, are clearly unresolved. Therefore, when we take into account the uncertainty in the measurements for these targets, we can only determine an upper limit of the N band envelope extension.

All targets have visibilities larger than 0.8. This corresponds, for a 130m baseline at 8µm to extensions smaller than 6 mas assuming an uniform disk or 3 mas assuming a FWHM Gaussian. However, to estimate an accurate envelope extension, we need to include the contribution from the central star to the total flux (star+envelope) for each of systems studied. Table 3. Published effective temperature for our target stars.

A simple two component model

If we assume, for all our targets, that the N band emission originates from the central star and and the circumstellar envelope only, we can write :

V tot = V e F e + V ⋆ (1 -F e ), (1)
where V tot , V ⋆ , and V e , are the total, stellar, and envelope visibilities, respectively, and F e , is the ratio of the envelope flux to the total object flux.

We can estimate the stellar angular diameter using the spectral class of each object and their distance derived from Hipparcos parallax measurements (Hipparcos, the New Reduction van Leeuwen, 2007). As shown in Table 2, all targets have stellar diameters smaller than 1 mas. Thus, V ⋆ > 0.99 at 8µm and we can assume V ⋆ =1 in the entire N band for these stars.

Using this assumption, we can rewrite Eq.( 1) :

V e = V tot -1 + F e F e . (2)
Consequently, in order to calculate the envelope visibility V e we only need the relative envelope flux F e .

Reconstructing the spectral energy distribution

We chose to constrain F e by fitting the objects' spectral energy distribution (SED) using stellar models. For each target, we reconstructed its SED using data available from the VIZIER 1 database (see (1976) and [START_REF] Thompson | Catalogue of stellar ultraviolet fluxes van Leeuwen[END_REF], magnitudes are from Morel & Magnenat (1978), Monet et al. (2003), [START_REF] Zacharias | NOMAD catalog[END_REF], and 2MASS [START_REF] Cutri | 2MASS all Sky Catalog of point sources[END_REF], and IRAS measurements are used for the mid to far infrared. We note that Be stars are observed to exhibit photometric variations with timescales from a few months to several years with amplitudes up to 0.5 Mag (Harmanec, 1983) so that the SEDs can slightly vary with time. However, contemporaneous UV to mid-IR photometry are not available to compare with our VLTI/MIDI measurements so we use measurements from 0.2 to 100µm obtained at other times. The resulting reconstructed SEDs for our targets are plotted in Fig 2.

The stellar parameters and distance

In order to constrain the stellar contribution to the flux we need T eff , R ⋆ , and g eff for the central star. However, these stellar parameters for rapidly rotating Be stars are quite uncertain. There are several methods in the literature that are used to determine these stellar parameters. For example, Dachs et al. (1989) derived stellar parameters and extinctions of 46 Be stars from UBV 1 http://vizier.u-strasbg.fr/viz-bin/VizieR Fig. 2. Reconstructed SEDs for the 7 Be star systems observed with the VLTI/MIDI instrument are shown with gray solid lines. Diamonds represent measurements taken from the VIZIER database. The portions of the SEDs up to 1µm (solid black lines) are fitted using reddened Kurucz models with the stellar parameters given in Table 4. Dashed lines represent the wavelengths for which the envelope flux is calculated, i.e., at 8µm, and 12µm. and JHKLMN photometric measurements. Later, Dachs et al. (1990) determined the stellar parameters of 26 objects from an analysis of the Balmer decrement. Chauville et al. (2001) and Levenhagen & Leister (2006) use spectroscopic measurements of photospheric lines to calculate these parameters. Table 3 summarizes T eff derived by all of these authors for our target stars.

The discrepancy between T eff determined with these methods for classical Be stars can reach several thousand of Kelvin for some targets (see ζ Tau for example). These ranges mainly stems from the latitudinal-dependent surface temperatures induced by rapid rotation. This effect will depend on the inclination angle of the star. Also, the circumstellar gas can cause reddening and obscure the star affecting the values of the spectrophotometrically determined parameters.

Fremat et al. ( 2005) accounted rapid rotation in their analyzes and quoted "apparent" stellar T eff and g eff , and then compared these values to non-rotating counter-parts that are typically thought to be ∼ 5 to 15% hotter. We adopt these apparent or average values for our work. The apparent T eff derived by Fremat et al. (2005) for our sample stars are also presented in Table 3.

Assigning a value for the stellar radius for a rapidly rotating star can also be difficult since the polar radius and equatorial radius can be quite different. For our targets, typical values for the stellar radius are deduced from the objects' spectral class and then adjusted to obtain the best SED fit in the UV and visible. The fit of the stellar flux is accomplished using Kurucz models [START_REF] Kurucz | [END_REF]. We note that the use of plan parallel atmospheric models with a global T eff for fast rotating star may not apply. However, the effect on the mid-infrared flux should be negligible for these wavelengths since the stellar radiation is dominated by black-body emission. The Kurucz models are reddened using the law of extinction from Cardelli et al.(1989) combined with the standard interstellar R v value of 3.1. Thus, with T eff and log g taken from Fremat et al. (2005), and the distance derived from Hipparcos parallaxes measurements, only the values for R ⋆ and A v are determined from the fit of the SED.

Estimating the errors on the envelope relative flux

To estimate the errors in the relative envelope flux F e at 8µm and 12µm, we have to take into account the uncertainty in T e f f adopted from Table 3, the possible interstellar and/or circumstellar reddening of the UV and visible stellar radiation, and the fact that a part of the visible emission can originate from the circumstellar envelope (up to ∼ 0.5 Mag).

ζ Tau, is probably the most challenging target to model since its MK type is not well constrained: the published luminosity class ranges from II to IV and its spectral type from B2 to B4. Moreover, at 0.2µm an absorption bump is clearly evident in its UV spectra. This is a indicates the presence of interstellar and/or circumstellar extinction of the stellar emission. In Fig 3 we plot the SED of ζ Tau with two "extreme" models, a reddened B2V for which up to 40% of the visible flux comes from the envelope, and a B4III unreddened star fully dominated in the visible by the stellar emission.

We obtain errors for F e on the order of 0.06 at 8µm and 0.04 at 12µm. This simple analysis shows that the calculated N band relative envelope flux is not strongly dependent on the spectral class and associated stellar parameters, but that the main uncertainty comes from the envelope contribution in the visible.

We applied this method to all our seven targets to determine the envelope relative flux F e and its uncertainty. The parameters adopted and the results of this SED model fitting are presented in Table 4. The stellar SEDs are also over-plotted in Fig 2.

, T eff =22000K, R ⋆ =5R ⊙ , A V =0.2, R V =3.1)
for which up to 40% of the visible flux comes from the envelope, and a B4III unreddened star dominated by the stellar emission in the visible (dotted line, T eff =16000K, and R ⋆ =8.5R ⊙ ).

N band envelope extension

We determine the N band envelope visibility, V e , using Eq. ( 2) and values of F e from Table 4. To determine the mean visibility value, V tot , at 8µm and 12µm, we use a linear interpolation of the measured visibility in the same wavelength range, assuming that the measurements at each wavelength are independent. This method allows us to slightly decrease the errors on the interpolated visibility. Finally, we determine the envelope extension using an equivalent uniform disk diameter. The results, assuming these simple models, are given in Table ??. For convenience, the K band measurements from VLTI/AMBER (Meilland et al. 2007a(Meilland et al. , 2007b(Meilland et al. , 2008) ) or CHARA (Gies et al. 2007) are also quoted if available.

As mentioned previously, the envelope is resolved only for α Ara, δ Cen, and ζ Tau (only at 12µm). For κ CMa, p Car, ζ Tau (at 8µm), α Col, and β CMi, we determine a upper limit to the envelope extension only. Please see the next section for more discussion about α Area's N and K band geometry, and on the compactness and wavelength dependence of the envelope size.

Discussion

The K and N band geometry of α Ara

First VLTI/MIDI data for α Ara were obtained by Chesneau et al. (2005). Unfortunately, the baselines used were too short and Gies et al. (2008 4 ) the data quality was not good enough to accurately constrain its geometry. The authors were only able to put an upper limit on the envelope extension, i.e. less than 10 mas at 8µm. However, this was already smaller than what Stee & Bittar (2001) predicted for this Be star. Thus, these authors proposed a possible disk truncation and found spectroscopic evidence for the presence of a companion orbiting outside the estimated outer-radius of the equatorial disk, confining the disk within the Roche lobe of the star.

Using VLTI/AMBER data, Meilland et al. (2007a) investigated α Ara K-band geometry and kinematics. They clearly resolved the equatorial disk and measured its flattening and orientation. The disk orientation was in agreement with the value deduced from polarimetric measurements by Yudin & Evans (1998) (P.A. =166 o ). Meilland et al. (2007a) also found that a few percent of the K band flux originates from an extended structure in the polar direction. Using the SIMECA code developed by Stee et al. (1996), they were able to model this object assuming a quasi-critical rotator surrounded by both a dense Keplerian rotating disk and an enhanced polar wind with terminal expansion velocities up to 1000 kms -1 . Note that this is the first time that the nature of the disk rotation, found to be Keplerian, was directly measured for a Be star.

Due to the SCI-PHOT mode and the use of longer baselines (i.e. UT1-UT4), the quality of our new VLTI/MIDI data was good enough to accurately constrain the envelope geometry in the N band. We were able to successfully model our measurements using a uniform ellipse, and found the same value for the major-axis as determined in the K-band from the VLTI/AMBER data by Meilland et al. (2007a). The envelope extension (a) and flattening (the ratio of a/b) are given in Table 6.

Both the VLTI AMBER and VLTI/MIDI measurements produced similar disk sizes over the wavelengths, 2.1µm, 8µm, and 12 µm. The envelope extension is plotted as a function of the baseline orientations in Fig 4 . A SIMECA image in the K band from Meilland et al. (2007a) is also plotted in this figure to allow a direct comparison with the K and N band measurements.

One of the most striking results from our modeling is the fact that the equatorial disk extensions remain relatively constant between 2.1µm and 12µm. We originally expected that longer wavelengths would probe cooler regions and thus larger disk volumes, which is clearly not the case. However, small variations on the order of ∼ 2σ are detected between 8µm and 12µm along the direction of the major-axis. This wavelength independent disk size may be evidence for disk truncation by a companion, in agreement with the suggestion of Meilland et al. (2007a). However, a free-free and free-bound radiation dominated envelope with a non-truncated isothermal equatorial disk may also λ(µm) major-axis (mas) flattening 2.1 7.3±2 2.7±0.9 8 5.8±0.8 2.4±1.0 12 8.2±1.2 2.6±1.2 Table 6. α Ara equatorial disk major-axis and flattening.

give the same result depending on thermal structure of the disk. For example, studies that have determined self-consistent disk temperature distributions based on the assumption of radiative equilibrium, find that, on average, the temperature distribution is nearly constant up to ∼ 100 R * and is approximatively half the T eff of the central star. However, depending on the equatorial density near the stellar surface, large variations in temperature can occur especially in dense disks. For example, see Millar & Marlborough (1998) for thermal structure of pure hydrogen disks, Jones et al. (2004) for the effects of iron line cooling on temperatures, or Sigut & Jones (2007) for models that include a realistic solar type chemical composition. Thus, it may be difficult to propose a general scheme for each of our program Be stars since the disk density may vary significantly from star to star.

Another interesting result can be seen from 1. The VLTI/MIDI angular resolution is 4-6 times smaller than the VLTI/AMBER. The polar wind extends out to ∼ 6 mas and is almost fully resolved by VLTI/AMBER with 80m baselines, whereas, the VLTI/MIDI visibilities with 120m baselines are ∼ 0.8. Thus, the overall effect on the total visibility for the VLTI/MIDI instrument is reduced by a factor of ∼5.

2. If the poles are hotter than the equatorial regions because of the effects of gravity darkening, the relative contribution to the total flux of the polar wind can significantly decrease between 2.1 and 8µm. In order to estimate this effect, a full modeling of the star+disk+wind would be required which is well beyond the scope of this paper.

The geometry of the remaining program Be stars

As mentioned above, the envelopes of all observed targets are much smaller than predicted by Stee & Bittar (2001). In fact, none of our disks are larger than 10 mas. Using the angular diameters from Table 3 and parameters from 1, it follows that the baseline orientations of α Col, β CMi, ζ Tau, and p Car correspond to an intermediate angle between their putative major and minor axis whereas δ Cen and κ CMa are certainly observed along their major axis.

Three stars in our sample were also observed with VLTI/AMBER : α Ara, κ CMa, δ Cen and their correspond-ing extensions in the K band continuum are 7.3±2, <3.6, and 1.6±0.4 mas, respectively. ζ Tau was observed with the CHARA interferometer in the K ′ band (Gies et al. 2007) and modeled with an elliptical Gaussian with a FWHM extension of 1.79±0.07 mas.

For κ CMa, the N band measurements are not accurate enough to determine if its envelope size varies with wavelength. On the other hand, δ Cen and ζ Tau seem to show some size variations between the K and N bands. For ζ Tau, Carciofi et al. (2009) have shown in a recent paper that a density wave seems to be orbiting within the disk. Thus, the apparent extension of the disk may be different at different times.As previously discussed, this is contrary to the findings for α Ara as we showed that its circumstellar disk extension is relatively constant over the infrared wavelengths considered.

Overall, for the Be star systems considered here, the variations in disk size between the K and N bands seem relatively small. This is quite different from observations in these bands for the dusty disks of supergiant B[e] stars. For example, Domiciano de Souza et al. (2007) find that the interferometric observations of the circumstellar environment for the B[e] star, CPD-57 o 2874 in the K and N band are well fitted by an elliptical Gaussian model with FWHM diameters that varies linearly with wavelength. However, the N-band emission mechanisms are different for these two classes of stars, i.e. free-free emission for classical Be stars and optically thick thermal dust emission for B[e] stars. Nevertheless, comparing K and N band emissions of circumstellar envelopes of hot stars may provide clues, allowing the physics of these systems to be deciphered.

Conclusion

This first mid-infrared interferometric statistical study of classical Be stars shows evidence for the compactness of their circumstellar environment in the N band, smaller than 10 mas for the all Be stars in our sample. Interestingly, Gies et al. (2007) find that Be stars in binary systems generally have smaller disks. Further, Be star disk extensions do not seem to depend strongly on the observed wavelength between 8 µm and 12 µm. For α Ara, an even stronger case can be made since the disk exhibits the same extension and flattening in the K and N bands. A possible result, that needs further investigation, is whether or not the disk size in various bands correlates with the T eff of the central star. Our investigation seems to suggest that there may be such a correlation, but this must be taken with caution considering the large uncertainties in the data and the fact that, in most cases, the baselines were not aligned with the disks' major axis.

This study can only be marginally extended to a larger number of targets because of VLTI/MIDI limited sensibility and accuracy. Fortunately, the number of targets available and the precision will significantly increase with MATISSE (Lopez et al. 2008), the second generation VTLI mid-infrared beamcombiner.

  Fig 2): UV measurements are from Jamar et al.

Fig. 3 .

 3 Fig. 3.A comparison of the reconstructed SED for ζ Tau (crosses) with two "extreme" models: a reddened B2V (solid line, T eff =22000K, R ⋆ =5R ⊙ , A V =0.2, R V =3.1) for which up to 40% of the visible flux comes from the envelope, and a B4III unreddened star dominated by the stellar emission in the visible (dotted line, T eff =16000K, and R ⋆ =8.5R ⊙ ).

  Fig 4. An apparent lack of an extended polar contribution in the N band data compared to the K band measurement is apparent. There are two possible explanations:

Fig. 4 .

 4 Fig. 4. α Ara envelope at 8µm (green bars) and 12µm (red bars) assuming a uniform disk + an unresolved star model plotted as a function of baseline orientation. VLTI/AMBER measurements (black bars) are from Meilland et al. (2007a) and the corresponding best SIMECA model in the continuum ( shaded background ) are also shown. The black, green, and red ellipses correspond to the uniform ellipse fit of the equatorial disk presented in Table 6 at 2.1, 8, and 12µm, respectively. The orientation of the major axis (dotted black line) is perpendicular to polarization angle measurement obtained by Yudin & Evans (1998).

Fig. 5 .

 5 Fig. 5. 8µm (black) and 12µm (gray) envelope extension as a function of T eff for the program Be stars. The dotted lines represent linear interpolations, accounting for the uncertainty of the measurements, of the temperature dependence on the N band size.

Table 1 .

 1 VLTI/MIDI observations for our program stars were carried out at Paranal Observatory between June 2004 and January 2009 with the 8.2 m UT telescopes. We obtained one visibility measurement for p Car, ζ Tau, κ CMa, and α Col, two for δ Cen and β CMi, and three for α Ara. Observations of 12 calibrator stars were also obtained during this campaign: HD40808, HD49161, HD50778, HD94452, HD 95272, VLTI/MIDI observing log for our program stars and their corresponding calibrators.

	Star Obs. Date Time	Baseline	Calibrators
	(UT) B n L(m) P.A.( o )	(HD)
	κ CMa 2008/12/14 02:33	117.0 25.2	40808, 94452
	δ Cen 2004/06/06 00:43 B 2 91.2 42.5	107446
	2008/12/13 08:27 B 1 129.2 27.7	40808, 94452
	α Ara 2006/04/15 03:50 B 1 126.8 1.7 95272, 129456, 151249
	2006/04/16 07:49 B 2 129.6 54.4	151249, 163376
	2006/07/08 04:15 B 3 120.7 74.1 169916, 198048, 217902
			218594
	p Car 2008/12/13 07:44	129.5 40.8	40808, 94452
	ζ Tau 2009/01/11 03:36	127.3 67.9	49161, 50778
	α Col 2008/12/13 05:17	128.5 64.0	40808, 94452
	β CMi 2008/12/17 07:07 B 1 129.5 63.8	49161, 50778
	2009/01/11 04:44 B 2 125.3 64.5	49161, 50778

Table 2 .

 2 Distance, spectral class, and estimated stellar angular diameters (Φ ⋆ ) for the seven observed Be stars.

	Name Distance (pc) Spectral Class κ CMa 202 ± 5 B1.5IV δ Cen 127 ± 8 B2IV	Φ ⋆ (mas) 0.44 ± 0.15 0.65 ± 0.3	Fig. 1. VLTI/MIDI visibilities of the 7 Be stars plotted as a func-tion of wavelength across the N band. For stars with multiple
	α Ara	82 ± 6	B3V	0.88 ± 0.2	measurements (i.e., α Ara β CMi, and δ Cen) the correspond-
	p Car	148 ± 9	B4V	0.27± 0.15	ing baselines B 1 , B 2 , and B 3 are plotted in black, gray, and light
	ζ Tau	136 ± 15	B4III	0.65 ± 0.2	gray, respectively. See Table 1 for description of the baselines.
	α Col	80 ± 2	B7IV	0.58 ± 0.15
	β CMi	50 ± 1	B8V	0.56 ± 0.1

Table 4 .

 4 Stellar parameters adopted in this work and resulting 8 and 12µm envelope flux.

		κ CMa	δ Cen			α Ara	p Car	ζ Tau	α Col	β CMi
			B 1	B 2	B 1	B 2	B 3		B 1	B 2
	φ (mas) at 8 µm	<5.2	4.9 ±1.8 4.0 ±1.8 2.7 ±1.1 3.9 ±1.4 5.5 ±0.3 <6.0	<2.6	<2.3	<6.8	<1.1
	φ (mas) at 12µm	<8.4	5.3 ±4.3 6.9 ±2.7 3.6 ±1.8 6.0 ±2.3 8.1 ±0.6 <9.5	5.7±2.2	<3.6	<10.5 <1.8
	φ (mas) at 2.1µm	<3.6 2	1.6±0.4 3			7.3±2 1	-	1.79±0.07 4	-	-

Table 5 .

 5 Envelope size from our VLTI/MIDI data and literature :Meilland et al. (2007a 1 , 2007b 2 , 2008 3 ) and

  Table 2, each disk is plotted (in R ⋆ ) in the N band as a function of the central stars' T eff in Fig 5. This figure seems to indicate that the envelopes of late type Be stars (α Col and β CMi) are smaller than those of earlier types. However, this claim should be taken with caution considering the uncertainty of the measurements.Moreover, with the exception of α Ara, one has to keep in mind that these stars are measured with baselines not necessarily aligned with their major-axis. Comparing polarization measurements of Yudin & Evans (1998) (see Table7) and the baselines orientations from Table

Table 7 .

 7 Yudin & Evans (1998) polarization measurements and baselines orientation with respect to the disk putative major-axis.

		κ CMa	δ Cen		α Ara		p Car ζ Tau α Col	β CMi
	Polarization(%)	0.3	0.33		0.58		1.16	1.13	0.15	0.14
	Polarization angle ( o )	106	137			166		68	34	109	90
	Baseline orientation ( o )		B 1	B 2	B 1	B 2	B 3		B 1	B 2
	compared to the disk major-axis	9	4.5 19.3 59 21.6 1.9	62.8	56.1	45	63.8 64.5

Acknowledgements. The Programme National de Physique Stellaire (PNPS) and the Institut National en Sciences de l'Univers (INSU) are acknowledged for their financial support. A. Meilland acknowledges financial support from the Max Planck Institut fur Radioastronomy, and thanks Thomas Driebe for his help in reducing MIDI data. CEJ wishes to acknowledge financial support from NSERC, Natural Sciences and Engineering Research Council of Canada. The authors would like to thank the referee, Thomas Rivinius, for his careful reading and suggestions that greatly improved the paper. This research has made use of SIMBAD and VIZIER databases, operated at CDS, Strasbourg, France.


Camille Guinaudeau

Andreu Girbau

Andreu Girbau Xalabarder

Textual Analysis for Video Memorability Prediction

Introduction

Efficiently predicting the memorability of a video could be very useful in the context of multimedia document analysis, e.g. automatic summarization. To progress on this multimodal task, the mediaeval evaluation benchmark organizes for the fifth year the Predicting Video Memorability task [START_REF] Sweeney | Overview of the MediaEval 2022 predicting video memorability task[END_REF]. Previous work on video memorability (VM) showed that semantic information plays a major role for VM prediction, either through manual descriptions ( [START_REF] Cohendet | Annotating, understanding, and predicting long-term video memorability[END_REF][START_REF] Cohendet | VideoMem: Constructing, analyzing, predicting short-term and long-term video memorability[END_REF][START_REF] Kleinlein | Topic-Oriented Text Features Can Match Visual Deep Models of Video Memorability[END_REF]), manual or automatic captioning ( [START_REF] Newman | Multimodal memorability: Modeling effects of semantics and decay on video memorability[END_REF], [START_REF] Cohendet | VideoMem: Constructing, analyzing, predicting short-term and long-term video memorability[END_REF], [START_REF] Shekhar | Show and recall: Learning what makes videos memorable[END_REF]) or semantic concepts [START_REF] Wang | RUC at MediaEval 2019: Video Memorability Prediction Based on Visual Textual and Concept Related Features[END_REF].

While they all agree that semantic information is important in distinguishing memorable videos from less memorable ones, they struggle to identify the semantic cue that helps make this distinction. [START_REF] Gupta | Linear Models for Video Memorability Prediction Using Visual and Semantic Features[END_REF] analyzed the vocabulary terms corresponding to the most positive and negative correlation coefficients of their caption-based models and showed that the highest negative coefficients are all dominated by terms related to natural scenery; whereas the highest positive coefficients are dominated by terms related to people. Similarly, [START_REF] Kleinlein | Topic-Oriented Text Features Can Match Visual Deep Models of Video Memorability[END_REF] showed that videos that are the closest related to nature and landscapes generally exhibit worse average memorability ratios, but found that in the Memory10k dataset most documents deal with people, and couldn't find a clear distinction between topics that are memorable and topics that are not. Finally, following the idea that people pay more attention to the concepts they are familiar with, [START_REF] Wang | RUC at MediaEval 2019: Video Memorability Prediction Based on Visual Textual and Concept Related Features[END_REF] used 156 familiar concepts [START_REF] Dolch | A basic sight vocabulary[END_REF] to generate ConceptNet feature vectors to represent the videos. However, the fusion of this semantic information failed to improve the results.

In this paper, we present an analysis of the multimodal characteristics that can be used to predict video memorability. First, we describe the characteristics and models used for the first subtask of the Prediction Video Memorability task on the Memento10k dataset [START_REF] Newman | Multimodal memorability: Modeling effects of semantics and decay on video memorability[END_REF]. We then analyze the results and study the videos' textual representations with respect to their memorability scores. Since manual textual representations may include a translation of whether humans find the videos memorable, we conduct these analyzes on both manual and automatic textual representations.

Video Memorability Prediction

Features To represent the videos, two different kinds of textual features were used. First, the manual descriptions provided by the organizers were represented by a 384 dimensional dense vector using Sentence-BERT [START_REF] Reimers | Sentence-BERT: Sentence Embeddings using Siamese BERT-Networks[END_REF]. Second, two automatic captions for each video were computed using the ClipCap Image captioning model [START_REF] Mokady | ClipCap: CLIP Prefix for Image Captioning[END_REF]. To do so, one keyframe was automatically extracted for each video using the Katna tool 1 , and the ClipCap model was applied using the model trained either on the MS-COCO dataset [START_REF] Lin | Microsoft coco: Common objects in context[END_REF] or on the Conceptual Caption dataset [START_REF] Sharma | Conceptual captions: A cleaned, hypernymed, image alt-text dataset for automatic image captioning[END_REF]. Captions from both models were concatenated and, similarly to the manual descriptions, were represented by 384 dimensional dense vectors using Sentence-BERT. In addition to these textual features, two visual features provided by the organizers were also employed: Densenet [START_REF] Huang | Densely connected convolutional networks[END_REF] and Resnet [START_REF] He | Deep residual learning for image recognition[END_REF].

Models and scores fusion

Two models were defined for memorability scores prediction. The first one, used for prediction based on textual representations, is a 3-layer sequential model. In order to avoid overfitting a dropout regularization is added on the visible layer. We use Huber loss as the loss function and a learning rate of 1 × 10 -5 . The second model uses video features extracted from deep learning models provided by [START_REF] Sweeney | Overview of the MediaEval 2022 predicting video memorability task[END_REF], over 3 frames -first, middle, and last-. We concatenate the features of the 3 video frames and input the resulting feature vector to a 3-layer sequential model. We use dropout on each layer, as loss function the Mean Squared Error (MSE), and a learning rate of 2 × 10 -3 . To keep the output values scaled from 0 to 1, we use a sigmoid function in the output layer for both models. Also, for both models, we make use of the Adam optimizer [START_REF] Kingma | A Method for Stochastic Optimization[END_REF]. To combine textual and visual features for memorability scores prediction, the embeddings are first fed to the models and the fusion scores 𝑆 𝐹 are then computed thanks to 𝑆 𝐹 (𝑖) = 𝛼 * 𝑆 𝑀 (𝑖) + (1 -𝛼) * 𝑆 𝑁 (𝑖), with 𝛼 ∈ [0, 1] and where 𝑆 𝑀 (𝑖) is the score for the 𝑖 𝑡ℎ video according to modality 𝑀 and 𝑆 𝑁 (𝑖) the score for the 𝑖 𝑡ℎ video according to modality 𝑁 .

Results and Analysis

Table 1 summarizes the results obtained on development and test sets of the Memento10k dataset. Concerning the scores obtained from textual representations, although the manual descriptions provide better scores, automatic captions give encouraging results, meaning that automatic captions carry information, automatically extracted from images, related to memorability. Combination of both textual features further improves the results, which suggests that the information carried by manual descriptions and automatic captions is not strictly redundant.

Visual characteristics provide correlation scores almost on par with those obtained from manual descriptions, and the combination of the two visual features (DenseNet and ResNet) increases the results. As before, the combination of visual features and automatic captions further improves the correlation scores. Finally, from this Table, we can see that, although the combination of automatically calculated representations (captions and visual features) provides encouraging results (0.62), the best results are still obtained when manual descriptions are used, suggesting that humans have the ability to introduce information related to memorability into their representations, that is not carried by visual characteristics or automatic captions.

Once again, these experiments show that semantic information is a mainstay in the prediction of video memorability. In order to understand what makes video memorable, we have analyzed the textual information attached with each video. First, we dig into topic analysis to verify if some concepts are more memorable than others, as claimed by [START_REF] Wang | RUC at MediaEval 2019: Video Memorability Prediction Based on Visual Textual and Concept Related Features[END_REF]. To do so, we identified the 10 words in the manual descriptions that have the most important differences in normalized frequency between videos associated with the higher memorability scores and the lower memorability scores. We then computed the normalized frequency of this 10 words for 3 splits of the development and training set regarding the values of the memorability scores: between 0.4 and 0.62 , between 0.6 and 0.8 and between 0.8 and 1.

Figure 1 shows the values of these normalized frequencies. Similarly to [START_REF] Kleinlein | Topic-Oriented Text Features Can Match Visual Deep Models of Video Memorability[END_REF], we couldn't see a clear separation in the topics addressed by the memorable videos and by the other videos. If, as presented in [START_REF] Gupta | Linear Models for Video Memorability Prediction Using Visual and Semantic Features[END_REF], the vocabulary associated with nature is more frequent in videos associated with a lower similarity scores, the human topic is addressed by all the videos, whatever their score. However, the way this topic is described in the manual descriptions is different between memorable videos (where the vocabulary is more precise man, woman, baby, girl) and non memorable ones, where the term used is more vague (people).

A possible explanation for this phenomenon is that annotators, remembering the video better, use more precise vocabulary in video with high memorability scores when writing the manual description. Another possibility is that the video is more memorable because the elements within it are more easily identifiable by viewers. This hypothesis is in line with the notion of familiarity mentioned by [START_REF] Wang | RUC at MediaEval 2019: Video Memorability Prediction Based on Visual Textual and Concept Related Features[END_REF]. To confirm this hypothesis, we carried out textual analyzestrough part-of-speech tagging and dependency parsing using the nltk tool [START_REF] Bird | Natural language processing with Python: analyzing text with the natural language toolkit[END_REF] -of the manual descriptions and the automatic captions, to check whether the sentences used to describe the most memorable videos are more precise than those used for the less memorable ones.

First, we found that the length of both manual descriptions and automatic captions were longer for memorable videos than for non-memorable videos, suggesting that either 1) memorable videos have more to describe or 2) elements present in the video are described more extensively for memorable videos. Moreover, we observed, as showed in the left part of Figure 2, that the proportion of nouns used in the texts (descriptions or captions) associated with the videos decreases when the memorability score increases while the proportion of adjectives grows. This observation leads us to believe that textual representations tend to be more accurate for the most memorable videos. To corroborate this assumption, the proportion of 8 frequent adjectives for physical descriptions and 8 common English colors3 were computed on manual descriptions and automatic captions. The right side of Figure 2 shows the evolution of this proportion and confirms that the proportion of colors and physical adjectives usage improves with the values of memorability scores on both manual descriptions and automatic captions. Finally, similar analyzes were conducted on the VideoMem dataset and the automatic captions associated with the videos follow the same pattern than those of Memento10k: length and adjectives proportion growing and nouns proportion decreasing while the memorability scores rise.

Future work

In this paper, we analyze the results of two sequential models based on visual and textual features to understand what makes a video memorable. We show that automatic captions and visual features can provide encouraging results but still miss memorability related information carried by manual descriptions. We also conduct an analysis of the manual and automatic textual representation of videos, showing that more memorable videos are associated with more precise descriptions even in automatic captions. In order to confirm this finding, it would be interesting to estimate the vagueness of the descriptions associated with videos, with a tool such as VAGO [START_REF] Guélorget | Combining vagueness detection with deep learning to identify fake news[END_REF], and use this vagueness prediction as a clue for memorability prediction.

Figure 1 :

 1 Figure 1: Normalized frequency of words regarding memorability scores.

( a )

 a Nouns and adjectives (b) Color and physical adjectives

Figure 2 :

 2 Figure 2: Proportion of adjectives and nouns in video descriptions w.r.t. memorability scores.

Table 1

 1 Spearman (Sp.) and Pearson (P.) correlation and Mean Square Error (MSE) values on development and test sets. The parameters column gives the weight values for the linear combination of scores, optimized on the development set.

	Runs	Param.		Development		Test
			Sp.	P.	MSE	Sp.	P.	MSE
	Desc.	-	0.599 0.599 0.008 0.571 0.575 0.007
	Captions	-	0.511 0.501 0.008	-	-	-
	ResNet (RN)	-	0.557 0.550 0.008	-	-	-
	DenseNet (DN)	-	0.561 0.550 0.008	-	-	-
	RN + DN	0.5	0.587 0.576 0.008 0.546 0.531 0.008
	Captions + (RN + DN)	0.3	0.622 0.610 0.007	-	-	-
	Desc. + Captions	0.6	0.629 0.620 0.007 0.589 0.594 0.007
	Desc. + (RN + DN)	0.5	0.668 0.666 0.007 0.626 0.628 0.006
	Desc. + Captions + (RN + DN) 0.7, 0.6 0.674 0.667 0.007 0.629 0.633 0.006

https://katna.readthedocs.io/en/stable/index.html

There is almost no video with scores lower than 0.4.

small, large, long, tall, little, big, young and old / white, yellow, green, blue, purple, red, orange and black.


Emna Soyed
email: emna1.soyed@gmail.com

Ramzi Chaieb
email: ramzi.chaieb@hotmail.com

Karim Kalti
email: karim.kalti@gmail.com

Superpixel Based Segmentation of Historical Document Images Using a Multiscale Texture Analysis

Keywords: Segmentation of Historical Document Images, Multiscale texture analysis, SLIC superpixel, Gabor features, Merging classification results

In this paper, a superpixel based segmentation of Historical Document Images (HDIs) using multiscale texture analysis is proposed. A Simple Linear Iterative Clustering (SLIC) superpixel technique and Kmeans classifier are applied in order to separate the input image into background and foreground superpixels. The foreground superpixels are characterized by the standard deviation and the mean of the Gabor features. These features are extracted in a multiscale fashion to adapt to the variability of the textures that may be present in HDIs. Text/graphic separation is then performed by applying a classification of the foreground superpixels for each texture analysis scale followed by a merging step of the obtained classification results. Since the classification results depend on the used classifier, a comparative study is performed for supervised (Support Vector Machine (SVM), K-Nearest Neighbors (KNN)) and unsupervised (Kmeans, Fuzzy C-Means (FCM)) techniques. Experiments show the effectiveness of our proposed method especially when compared with similar work in the literature.

Introduction

In recent years, digitizing collections of cultural images becomes the vision of libraries and museums to ensure a sustainable conservation of historical collections and a world access to large sets of cultural heritage documents. Thus, many open issues and challenges have been raised, such as developing content-based image indexing and retrieval tools and ensuring the efficiency of document interpretation systems [START_REF] Loussaief | Machine learning framework for image classification[END_REF][START_REF] El Bazzi | A graph based method for Arabic document indexing[END_REF]. In particular, the literature shows that a huge amount of issues is related to the particularities of historical documents such as noise, degradation, page skew, specific fonts, irregular spacing between characters, random alignment and different text orientations [START_REF] Coustaty | Historical document analysis: a review of French projects and open issues[END_REF]. Recently, the family of superpixel based algorithms has gained a great attention of many researchers in how we can use segmentation methods to gain new insights into the world of image processing. The superpixel technique uses the degree of characteristics similarity between pixels and acquires the redundant information in the image to grip pixels. Therefore, they can reduce the complexity of image post processing tasks and by following will lead to a more effective treatment.

Related Work

For instance, Jian used the superpixel technique to perform the image background segmentation [START_REF] Jiang | Linear solution to scale invariant global figure ground separation[END_REF]. The particular aggregation of information provided by the superpixel has been proved by Li et al. to be useful for image segmentation [START_REF] Li | Segmentation using superpixels: a bipartite graph partitioning approach[END_REF]. Cohen et al. used spatial and color features extracted from superpixels to separate drawings regions from background of ancient documents [START_REF] Cohen | Robust text and drawing segmentation algorithm for historical documents[END_REF]. Garz et al. used a multistage algorithm based on interest points to detect the layout entities in ancient manuscripts [START_REF] Garz | Layout analysis for historical manuscripts using sift features[END_REF]. They proposed a new stroke width computing method by using intensity and Simple Linear Iterative Clustering (SLIC) superpixels region growing to segment text from low quality images. The limitation of this method is the fact that it cannot be used in text images with various strokes and it was proved to not be good for complex background text segmentation. Mehri et al. pointed out a background separation technique for ancient documents using spatial and color features extracted from superpixels [START_REF] Mehri | Use of SLIC superpixels for ancient document image enhancement and segmentation[END_REF]. Moreover, in order to extract the text and graphic regions, they presented an algorithm based on the use of the SLIC superpixel approach and Gabor descriptors to segment historical document images [START_REF] Mehri | Learning texture features for enhancement and segmentation of historical document images[END_REF]. Mehri et al. used different sliding window sizes for the texture analysis. The descriptors extracted from these different windows were embedded in a single vector which may affect the performance of the classification step.

In this paper, we propose a contribution to overcome these inefficiencies by developing a HDIs segmentation method based on using multiscale texture analysis. The proposed method is based on the SLIC superpixel technique and Gabor features. Accordingly, four different scales analysis are used and four classifiers are applied to extract textual information from the graphical ones. In the literature, many methods for text/graphic classification were elaborated [START_REF] Mehri | Use of SLIC superpixels for ancient document image enhancement and segmentation[END_REF][START_REF] Mehri | Learning texture features for enhancement and segmentation of historical document images[END_REF]. Unfortunately, all these works did not make any comparative study between them. Our objective is to apply a comparison between four different classification methods in order to ensure the effectiveness of the selective classifier. First, an interactive feature learning step is introduced to train supervised classifiers. Then, a merging method is developed to deal with the variety of the multiscale analysis windows and to make easier the choice of the best classifier.

The remainder of this paper is organized as follows. Section 3 presents an overview of our proposed method. The adopted superpixel technique is described and the investigated Gabor features are detailed. In order to outline the pertinence of the experimental protocol, Sect. 4 describes the ground truth, the segmentation results and their evaluation. Finally, Sect. 5 gives a conclusion in the end of this paper.

The Proposed Method

In this work, we are not searching for an accurate segmentation, but to find regions with similar characteristics in a variety of scales and fonts, and to discriminate between text and graphic components. Figure 1 presents the schematic block representation of the proposed method.

The proposed method starts by separating the input HDI into foreground and background superpixels using the SLIC superpixel technique and the Kmeans clustering to obtain an enhanced and non noisy background and to avoid the complexity of pixel based segmentation. The extracted Gabor features for each texture analysis scale are fed as input to four different classifiers in order to classify the foreground superpixels into two classes (text or graphic). Finally, a merging step is introduced in order to reduce the multi-variance of the sliding windows on the classification step. In order to have a relevant segmentation of textual regions from graphical ones and a separation between different Kinds of graphics and various text fonts, the proposed method is composed of four steps: Preprocessing, Texture feature extraction, Superpixels classification and Merging results.

Preprocessing: Foreground/Background Segmentation

First, a HDI is fed as input and converted to a grayscale image. The proposed method does not assume a prior knowledge regarding document content and layout. In order to improve the segmentation results, a median filter is applied on the input image. Figure 2 shows an example of the preprocessing step. Once the preprocessing is terminated, the image is segmented. The aim of the segmentation task is to simplify the representation of an image on some things easier and more meaningful for the extraction and analysis step. Rather than using a rigid structure of pixel grid, the superpixel technique has been used to group pixels sharing similar properties or characteristics into a significant polygon shaped region. The superpixels are used in our method as basic units instead of using the rigid structure of pixel grid and become a consistent alternative for a foreground/background segmentation. Achanta et al. performed an empirical comparison of five state-of-the-art superpixel methods [START_REF] Achanta | SLIC superpixels compared to state-of-the-art superpixel methods[END_REF]. They classified the existing superpixel methods into three classes: gradient-based-ascent, graph-based and SLIC superpixel. They stated that the SLIC technique is the best superpixel method regarding the segmentation results and performances, the ability to adhere to image boundaries, the efficiency of memory, the simplicity of use and the ability to control superpixel compactness. Actually, by setting the superpixels number equal to 0.08%, the SLIC superpixel technique is performed on a grayscale image. The segmented image is generated by grouping pixels sharing gray level similarity and their proximity in the image plane into significant polygon-shaped regions. Afterwards, the background and foreground clustering task is performed based on the calculation of the mean gray level value of each superpixel by averaging over all gray level pixels associated to the superpixel regions and using the Kmeans algorithm. By setting the number of clusters equal to two, two well separated clusters are extracted: one for the foreground and one represents the background information. An enhanced background is obtained by allocating the value of a white pixel to the superpixel center and the pixels belonging to them.

Texture Feature Extraction

The texture features explored in this work are based on a multiscale Gabor filter. The choice of Gabor features was based on an empirical comparison between texture features extracted using Gabor filter and the discrete Fourier transform presented by Ursani et al. [START_REF] Ursani | Texture features based on Fourier transform and Gabor filters: an empirical comparison[END_REF]. They proved that Gabor filters have good performance as opposed to Fourier transform, due to the optimal localization properties to get information in both frequency and spatial domains from the analyzed images. Moreover, some of the known texture-based approaches have been compared by Mehri et al. [START_REF] Mehri | Texture feature evaluation for segmentation of historical document images[END_REF]: The auto-correlation function, Gray level Co-occurrence Matrix, and Gabor features, used for the segmentation of digitized simplified ancient document images. They concluded that Gabor features results, obtained on ancient document images, are the best for distinguishing textual regions from graphical ones and for font segmentation. Nevertheless, few studies have been sought for ancient document image segmentation using the multiscale Gabor filters [START_REF] Raju | Text localization and extraction from complex color images[END_REF]. Therefore, the investigated Gabor features tested in our work are Gabor filters.

Gabor features are produced at several frequencies and orientations for texture characterization. Four orientations and six frequencies are used in common implementation [START_REF] Charrada | Texture approach for nets extraction application to old Arab newspapers images structuring[END_REF][START_REF] Zhong | Image patches analysis for text block identification[END_REF]. In this study, 24 Gabor filters are applied by using 4 different orientations {0, π 4 , π 2 , 3π 4 } and 6 distinct frequencies

{2 √ 2, 4 √ 2, 8 √ 2, 16 √ 2, 32 √ 2, 64 √ 2}.
The energies and the amplitudes of the Gabor filter outputs are investigated. Gabor features are extracted by the convolution of the analyzed document image and pairs of orientations and frequencies. The standard deviation and the mean values of Gabor filtered magnitudes and energies responses corresponding to all pixels defined in the same superpixel are extracted. The convolution of the gray level ancient document image with 24 Gabor filters are applied for each pixel defined in the same superpixel. After all, an easy way to extract Gabor features on the hole transformed image by the selective Gabor filter, is the use of a multiscale analysis technique to extract Gabor features from the selected foreground superpixels at four different sliding windows ((8 × 8), (16 × 16), (32 × 32) and (64 × 64)). Thus, four feature vectors are formed through four different sizes of sliding windows.

Superpixels Classification

Once the extracted Gabor features step is terminated, a classification task is performed by partitioning Gabor based features sets into separated classes. In order to select the best classifier for the discrimination of textual content from graphical one, a comparative study based on using several classifier techniques (Kmeans, Fuzzy C-Means (FCM), K-Nearest Neighbors (KNN), Support Vector Machine (SVM)) is carried out. Supervised and unsupervised techniques have been used on partitioning the image. As Kmeans and FCM are unsupervised techniques, only the number of classes K is provided as input. By setting the number of classes K equal to two, two well separated classes are extracted: one for textual region and one for graphical one. The SVM and KNN are supervised techniques used with a default kernel, in addition to the number of classes K, models are trained to generate the classification. Therefore, the training data is obtained by a manual selection of six regions (three belong to the text regions and three belong to the graphical ones) built for each image. Our goal is to construct a training data for each HDI, to ensure the effectiveness of the interactive method. The foreground superpixels issue is presented as a binary classification: one class for textual superpixels and another one containing all superpixels considered as other contents. Each superpixel belonging to the training data is noted as belonging to one of the two classes.

Merging Results

In many fields such as multiscale sliding windows, we are often confronted with multiple and conflicts sources of information. In this paper, we propose a merging method based on using the majority vote in order to resolve the conflicts of windows resolution. Four classification results are fed as input. Then, for each superpixel, we assign it to the majority voted class. If we have an equal case, we take the majority vote of neighboring superpixels already classified. Otherwise, we assign it to the same class of the nearest neighbor. A detailed schematic block representing the merging results of each Gabor filter application window for one classifier is illustrated in Fig. 3.

For each superpixel, the majority vote can be formed as follows:

V (s) = argmax i∈{1,2} (g(c i , s)) ( 1
)
Where s is the superpixel to be classified, V(s) is the majority vote and g(c i , s) is calculated as:

g(c i , s) = 4 j=1 w j ( 2
)
Where

w j = 1 if R j (s) = c i else w j = 0. R j (s)
is the classification result of s by the j classifier.

In this study, we consider only two classes which are text (c 1 ) and graphic (c 2 ). In the first equal case where g(c i , s) = 2, the majority vote among the neighboring superpixels is calculated based on the same definition. In the second equal case where g(c 1 , s) = g(c 2 , s), the superpixel can obtain its class from the nearest neighbor. For this purpose, the euclidean distance is used to select the nearest neighbor. This distance is calculated as follows:

dis s1,s2 = (x s2 -x s1 ) 2 + (y s2 -y s1 ) 2 (3)
Where s1 and s2 are two neighboring superpixels and (x s1 , y s1 ) and (x s2 , y s2 ) are respectively their coordinates of the centroids. Otherwise, we consider the result of the sliding window resolution 32 × 32 as the final decision class. The respective results of each step of our document image segmentation system will be reported and discussed in the next section.

Evaluation and Results

In this section, we perform experiments to illustrate and to discuss the performance of the proposed method. A brief description of the ground truth and the experimental corpus is presented.

Corpus and Ground Truth

To prove the influence of the merging step on extracting graphical regions from textual ones, we focus in our corpus on collecting HDIs containing various textual and graphical contents as shown in Fig. 4. A significant number of historical documents with both complex and simple layouts have been selected. This dataset is composed of grayscale, binary and color HDIs. Moreover, we have constructed the training data for each image by selecting three regions from the text and three others from the graphical ones. The training data is made for each image to guarantee a good understanding of the behavior of the extracted Gabor features.

Experiments and Results

The obtained results of the proposed segmentation method using the SLIC superpixel technique, multiscale Gabor filters and several classifiers techniques are illustrated in this section. By visual inspection of the obtained results, we note that the proposed method in this article ensure satisfying results particularly in distinguishing textual regions from the graphical ones. We perform a qualitative comparison of four different classification techniques. Arguably, the SVM classifier highlights the best classification results for separating textual regions from the graphical ones. Figure 5 shows an example of HDI segmentation assigned to each classifier. The quantitative evaluation of the proposed segmentation method is made by four measures: Complete (CM), Quality (CR), Accuracy (ACC) and F-measure (F) based on the Confusion Matrix decisions. Table 1 presents the obtained measures values computed for each classifier. Thus, by merging the Gabor features multiscale application results, we note that the SVM technique is more relevant than other techniques used in this work for the characterization and segmentation of textual regions and graphical ones. We note that the superpixels are classified with 94% (CM), 92% (CR) and 91% (ACC). Further, the SVM technique based on learning Gabor features is more efficient than the method presented in [START_REF] Mehri | Learning texture features for enhancement and segmentation of historical document images[END_REF] as shown in Table 1. Furthermore, it tends to well classify graphical superpixels into graphical classes thanks to the use of an interactive method to create the training data. As shown in Table 2, the obtained F-measure values are congruent and very promising.

Conclusion

In this paper, we proposed an interactive method based on using SLIC superpixels and Gabor features to extract textual and graphical regions in HDIs. In this study, each superpixel is characterized by the average and the standard deviation of the response of Gabor filters application. The features are computed for different sliding window sizes. Then, the extracted features are classified using supervised (KNN, SVM) and unsupervised (Kmeans, FCM) clustering techniques. In order to deal with the variety of windows, the results of the different classifiers are merged using the majority vote. Finally, the evaluation shows that our approach brings satisfactory results and outperforms the one proposed in [START_REF] Mehri | Learning texture features for enhancement and segmentation of historical document images[END_REF] thanks to the merging step. Our further work will be the use of the presented algorithm on an another type of database for example administrative documents. We will also focus on evaluating and merging other texture extraction methods.

Fig. 1 .

 1 Fig. 1. Schematic block of the proposed method based on texture features extraction for HDI and merging results of each classifier.

Fig. 2 .

 2 Fig. 2. Illustration of the results of the preprocessing step. (a) and (f) Original HDIs (b) and (g) Grayscale conversion results (c) and (h) Results of SLIC superpixel segmentation (d) and (i) Results of the foreground/background segmentation (foreground superpixels are labeled as blue and the background ones are labeled as green) (e) and (j) Foreground SLIC superpixels.

Fig. 3 .

 3 Fig. 3. Proposed merging method.

Fig. 4 .

 4 Fig. 4. Illustration of some examples of HDIs containing textual and graphical contents ((a), (b), (c) and (d)) or only textual contents ((e), (f), (g) and (h)).

-

  The Complete (CM): The CM measures the probability of having a pixel that belongs to both the system response and the Ground Truth.-The Quality (CR): The CR is described by the percentage of the region (text or graphic) correctly extracted. -The Accuracy (ACC): The ACC consists of the combination of the two other CM and CR measures. The ACC can be defined as the percentage of the text and graphic region extracted by the algorithm. -The F-measure (F): The F-measure is calculated as a score resulting from the combination of CM and CR. It assesses both the completeness and the homogeneity criteria of a clustering result.

Fig. 5 .

 5 Fig. 5. Example of resulting images of the proposed method assigned to each classifier: (a) Original HDI (b) The foreground/background separation results (c) The foreground superpixels (d) The ground truth generated by GEDI [16] (e) Kmeans results (f) FCM results (g) KNN results (h) SVM results.

Table 1 .

 1 Evaluation and comparison of the proposed method of HDI segmentation by calculating several measures FCM Kmeans SVM KNN Method proposed in[START_REF] Mehri | Learning texture features for enhancement and segmentation of historical document images[END_REF]

	CM 0.52 0.67	0.94 0.90 0.87
	CR 0.63 0.64	0.92 0.91 0.83
	ACC 0.47 0.49	0.91 0.88 0.79

Table 2 .

 2 Evaluation of HDI segmentation for the classifier SVM by computing the F-MEASURE values (F)

	Only one font One font and one graphic Total
	F 0.88	0.92	0.90


Lucas Gruss

Philippe Chevrel

Mohamed Yagoubi

Maxime Thieffry

Alain Grossetête

Moving horizon estimation of xenon in pressurized water nuclear reactors using variable-step integration

A specific moving horizon estimation scheme is proposed for reconstructing the axial xenon concentrations in pressurized water nuclear reactors (PWR). The model considered is continuous time, with a parcimonious spatial discretization. The estimation scheme involves continuous time integration while using low frequency sampling of measurements. The results support a sensitivity analysis. The quality of the estimation is shown to be good, even though the model is particularly stiff, and the implementation altogether light.

I. Introduction

A. Context

Operation of nuclear power plants (NPP) requires that many variables be kept within specific bounds in order to provide safe and efficient production of electricity. One such quantity is the build-up of xenon as a fission product. No measurement of xenon in the fuel is available.

Fission of 235 U produces mostly 135 I with rate γ I , which β --decays into xenon with half-life λ I . In smaller proportions uranium directly produces xenon with rate γ X . 135 Xe has a half-life of λ X = 9.2 hours and is also transmuted into stable 136 Xe upon absorbing a neutron with rate σ X , which makes it a neutron poison. The xenon/iodine equations are the following, as can be described in [START_REF] Reuss | Précis de neutronique[END_REF]:

İ(t, z) = γ I Σ f Φ(t, z) -λ I I(t, z) Ẋ(t, z) = γ X Σ f Φ(t, z) + λ I I(t, z) -(λ X + σ X Φ)X(t, z) . (1)
The terms Γ I = γ I Σ f and Γ X = γ X Σ f are introduced to simplify the expressions. Σ f depends on the cross-section of the fuel, which is dependent on the fuel cycle. Φ(t, z), X(t, z) and I(t, z) describe the densities in neutrons, xenon and iodine in the core as functions of time t and altitude z.

The interaction between xenon and neutron flux is to be considered in power plant operation in order to avoid a phenomenon called xenon oscillations. Such oscillations can be triggered by a small change in power distribution that alternatively causes xenon to accumulate in one part of the core and be depleted in another part. If left uncontrolled, such oscillations can diverge in a positive feedback loop and lead to safety issues.

Common ways to reduce the risks associated with xenon are appropriate control of the axial offset (AO) of power as well as sufficient downtime after shutdown. Generally, operators have to wait dozens of hours before powering up a plant again. Correct estimation of the xenon in the core can therefore provide valuable information to power producers for diagnostics, start-up and safety.

B. Contributions of the paper

In the presented work, we implement moving horizon estimation (MHE) in order to estimate the axial xenon concentration and distribution in a pressurized water reactor (PWR). The performance of MHE requires a continous-time system to be numerically integrated over "large" time intervals to account for frequently varying inputs and fast dynamics while simultaneously accomodating the slow xenon/iodine dynamics. We rely on a judicious choice of integrator to accelerate the numerical integration between the sampled measurements. This as opposed to a costlier evaluation of a high-frequency discretized model, and for better control of the integration error. We also extend previously used MHE implementations for the xenon estimation task, providing sensitivity analysis on the model parameters.

C. Related works

Estimation of xenon in nuclear reactors has been tackled both by academics and industry practitioners. The work in [START_REF] Ponçot | Variational assimilation for xenon dynamical forecasts in neutronic using advanced background error covariance matrix modelling[END_REF] is notable as it uses variational methods to estimate xenon in the core, which are essentially MHE methods developed by the meteorology community. It introduces two methods called "3D-Var" and "4D-var", which are moving horizon estimators with respective estimation horizon 1 and N > 1. The model uses 30 nodes core and includes neutronics, thermodynamics and xenon poisoning. In the same vein, we aim in the current paper to resume this work with the underlying objective of computation effort parsimony. Moreover, a sensitivity analysis of the estimator to parametric model uncertainty is proposed, which are not presented in the previous work.

The authors of [START_REF] Lin | Control of spatial xenon oscillations in pressurized water reactors via the kalman filter[END_REF] use Kalman filtering to carry joint state and parameter estimation on a two-points model of the core. In [START_REF] Ansarifar | Sliding mode observer design for a PWR to estimate the xenon concentration & delayed neutrons precursor density based on the two point nuclear reactor model[END_REF], a sliding mode observer is synthesized on the same kind of two-point model. Two-points models are the bare minimum in order to model AO of the core but are imprecise. Finer spatial discretization can help bridge the representativity gap, which is the point of the model we present. MHE proves to outperform Kalman filtering and sliding-mode observation, which motivates our choice of MHE over them.

In [START_REF] Park | Estimation of neutron flux and xenon distributions via observer-based control theory[END_REF] an infinite dimension Luenberger high-gain observer is designed. The authors use a linearized model of the onedimensional one-group diffusion equation and xenon/iodine dynamics about the steady-state xenon and flux distributions and manage to successfully reconstruct xenon, iodine and flux distributions. Our model is less conservative as we do not perform linearization. While it is interesting that the observer is designed on an infinite-dimensional model, Luenberger observers suffer from the same drawbacks as the Kalman filter and require rather strong assumptions to be applied on non-linear models.

The work in [START_REF] Yin | Distributed moving horizon state estimation of two-time-scale nonlinear systems[END_REF] addresses moving horizon estimation for singularly perturbed systems by taking a distributed approach, requiring the design of a MHE for the fast dynamics and a MHE for the "slow" dynamics. Our application is slightly different as we are not primarly interested in the fast dynamics (neutrons in our case) and focus instead on the slower dynamics of xenon. We implement and design a single estimator to estimate both slow and fast states, which is a less demanding design approach.

II. Model

We provide details on the model of the nuclear plant that will support the moving horizon estimation scheme, represented on figure 1. The reader can find more in-depth explanations on the model in [START_REF] Dupre | Design and comparison of two advanced core control sys-tems for flexible operation of pressurized water reactors[END_REF] and [START_REF] Lemazurier | Conception d'un système avancé de réacteur PWR flexible par les apports conjoints de l'ingénierie système et de l'automatique[END_REF] . The model was used for predictive control in the previous works and was validated to be a good compromise between representativity with respect to the true physics of the core and complexity. The model assumes one delayed neutron group and includes neutron diffusion equation, thermal equations, xenon/iodine dynamics and counter-reactions through moderator and doppler effects, control rods action, boron concentration and xenon poisoning. We include a modelisation of the chemical and volumetric control system (CVCS), which is responsible for the boration and dilution in the core. Fig. 1: Simplified model for the nuclear core, according to [START_REF] Dupre | Design and comparison of two advanced core control sys-tems for flexible operation of pressurized water reactors[END_REF] The model is discretized into 6 axial nodes. In each node a state comprising xenon concentration, iodine concentration, neutron power density and delayed neutrons is considered. We include temperatures between the nodes for a total of 7 temperatures, as well as the steam generator dynamics: cold and hot leg temperatures. Additionally, two boron concentrations are modeled : one for the CVCS system and one inside the core. The boron concentration in the core is considered to be uniform. The model has a total of 35 states.

Temperatures between meshes are obtained by assuming that all the power produced in the mesh is transfered to the coolant. This heat transfer is modeled through a firstorder model in each node. The neutron density obeys a diffusion equation, which is implemented by considering an exchange rate D between adjacent nodes. Delayed neutrons are produced as the result of the disintegration of the fuel and other fission products with half-life β. A single delayed group is considered. Model parameters have been calibrated with the help of the neutronics codes available at Framatome.

While reactivity is usually a macro-level term in the whole core, we introduce a reactivity term in each node, which impact the neutron density in each node. This reactivity term can be broken down into five sub-terms:

• doppler anti-reactivity: depending on the local power level, the probability of fission is impacted. • moderator anti-reactivity: a local increase in temperature in the coolant decreases the density of the water. Its ability to slow neutrons down impacts the probability of fission. • xenon anti-reactivity: local build-up of xenon in the fuel will result in poisoning of the node, which will consume a considerable amount of neutrons, thus decreasing the number of fissions. The parameter K X relates xenon concentration to anti-reactivity. • boron anti-reactivity: the diluted boron absorbs neutrons and limits fissions. • control rods anti-reactivity: the control rods contain neutron absorbing material. Inserted from the top, they locally limit the number of fissions that can happen. We have hinted on the stiffness of the model and will expand on this topic now. The model is an ordinary differential equation, to be simulated by numerical integration. Stiffness arises from the different time-scales in the physics of the system. Xenon dynamics is orders of magnitude slower than the neutronics in the core. Integrating this model requires to use a very short integration interval in order to simulate the system. We leverage variable-step integrators as a way to speed up implementation of the estimator and limit the number of terms in the cost function.

The measurements which will be used in the cost function are 6 axial neutron flux measurements n 1..6 from excore neutronics chambers, as well as the cold leg and hot leg temperatures (T out and T in ) of the steam generator, for a total of 8 measurements.

We consider the power extracted at the steam generator P turb , control rods position h rods and dilution and boration flowrates u D/B as inputs to the system.

III. Moving horizon estimation

A. Principles

Moving horizon estimation was proposed in [START_REF] Jazwinski | Limited memory optimal filtering[END_REF] back in 1968, making the case for least-square optimization to solve state estimation and improve robustness to model uncertainty compared with extended Kalman filtering [START_REF] Kalman | A new approach to linear filtering and prediction problems[END_REF]. Kalman filtering and its variants such as the extended and "sigma-point" filters [START_REF] Van Der Merwe | Sigma-point kalman filters for nonlinear estimation and sensor-fusion: Applications to integrated navigation[END_REF], [START_REF] Julier | New extension of the kalman filter to nonlinear systems[END_REF], [START_REF] Nørgaard | New developments in state estimation for nonlinear systems[END_REF] have been the go-to techniques for nonlinear estimation due to their relative ease of implementation and good practical performances. At the same time, the drawbacks are the difficulties to consider hard constraints on the estimates and in dealing with non-linearity. On the other hand, MHE can elegantly include these characteristics in the optimization problem that it solves online. Another practical advantage of MHE over other estimation techniques is explored in [START_REF] Dubois | Performance evaluation of a moving horizon estimator for multi-rate sensor fusion with time-delayed measurements[END_REF], in which MHE is shown to be an interesting choice for multi-rate and time-delayed systems.

An overview of MHE is available in [START_REF] Rawlings | Model predictive control: Theory, computation, and design[END_REF], this estimation method has seen increasing attention thanks to model predictive control (MPC) and its breakthrough among practitioners and academia. MHE is considered to be the dual of MPC and as such most of the progress on predictive control, such as software availability and execution speed has been carried over to MHE without extra costs. Consider the following non-linear state-space model and its discretization :

ẋ(t) = f (x(t), u(t)) y(t) = h(x(t), u(t)) ⇒ x k+1 = Φ(x k , u k ) y k = h(x k , u k ) . (2)
Moving horizon estimation is most often formulated in discrete-time so it is convenient to introduce a discretization of the model, which can be obtained from numerical or exact integration. At each time t k = k × dt, past measurements [y k-N+1 , ..., y k ] and inputs [u k-N+1 , ..., u k-1 ] are used to estimate all the past states in [t k-N+1 , t k ]. Hence, at time t k , N estimates of the state are produced. We introduce the notation x[k -i|k], ∀i ∈ [0, N -1] to designate the state estimates obtained from running the MHE at time t k and Xk = {x[k -N + 1|k] . . . x[k|k]} is the set of all estimates produced. x k-N+1 is a prior estimate of the beginning of the estimation window. The following optimization problem is solved online:

min Xk J = N-1 i=0 w i 2 Q + N i=0 v i 2 R + x[k -N + 1|k] -x k-N+1 2 S k st. w i = x[k -i + 1|k] -Φ( x[k -i|k], u k-i ) v i = y k-i -h( x[k -i|k], u k-i ) x[i] ∈ X, u[i] ∈ U, w[i] ∈ W, v[i] ∈ V . (3)
The estimation window is shifted forward in time as new measurements are available on the system, and a new optimization problem is solved. Three terms can be outlined from the cost function:

• The term N-1 i=0 w i 2
Q is optionnal depending on the implemented algorithm. It is used when state disturbances are considered. Confidence in the model is tuned through the matrix Q: if confidence is lower the contribution of this term in the cost function should be lower, and conversly.

• The second term N i=0 v i 2 R is used to compare the estimates with the available measurements on the system, weighted with matrix R which tunes the confidence in the measurements.

• The term x[k -N + 1|k] -x k-N+1 2
S k is called "arrival cost". Ideally, the data used for filtering would go consist in all past measurements available on the system. This estimation scheme is called Full Information Estimation (FIE) [START_REF] Rawlings | Model predictive control: Theory, computation, and design[END_REF]. In practice, both memory and computation requirements would grow too large to solve the problem in finite-time. Therefore in MHE, a finite-sized window of past measurements is introduced. The arrival cost can be seen as a way to "approximate" the rest of the cost function of FIE. Important results on stability and robustness of MHE depend on the design of the arrival cost, especially when the system is only detectable. We present a common implementation of the arrival cost. Consider x k-N+1 = x[k -N + 1|k -1], obtained at the previous MHE step as prior estimate. The term S k is given by running an extended Kalman filter in parallel, setting S k = P -1 k where P k is the solution to the Riccati Equation associated with the EKF. Since P k approximates the covariance cov( k ) of the estimation error k , S k gives the confidence level in the previous estimation. The usual limitations on the EKF apply here, since P k is only a local approximation of cov( k ) of the actual state. The other common issue of the EKF is that the estimation error is assumed to be gaussian whereas transformation of white noise by a nonlinear function is not gaussian in the general case. Other methods to compute the arrival cost are explored in [START_REF] Ungarala | Computing arrival cost parameters in moving horizon estimation using sampling based filters[END_REF], comparing EKF, UKF and particule filter. In [START_REF] Sanchez | Adaptive arrival cost update for improving moving horizon estimation performance[END_REF], the arrival cost is computed using adaptive estimation techniques. It should also be noted that stability of MHE can be ensured even with a choice of a constant value for S k (diagonal in the proof in [START_REF] Alessandri | Advances in moving horizon estimation for nonlinear systems[END_REF]).

Under mild assumptions on uncertainties and noise (boundedness and compactness) [START_REF] Alessandri | Advances in moving horizon estimation for nonlinear systems[END_REF], the estimation scheme shows good stability and robustness properties that make this method attractive despite the associated costs in computation. Efforts to reduce the computational burden leverage sensitivity results from non-linear programming, leading to implementations using suboptimal optimization [START_REF] Alessandri | Advances in moving horizon estimation for nonlinear systems[END_REF], realtime iteration [START_REF] Vukov | Real-time nonlinear MPC and MHE for a large-scale mechatronic application[END_REF] or "one-step" MHE [START_REF] Kã | A real-time algorithm for moving horizon state and parameter estimation[END_REF].

B. A formulation for stiff systems

The dynamics of xenon is slow, due to the half-lives of 135 Xe and 135 I in the range of several hours. Such time constants call for long estimation windows. In doing so, we can have significant contribution of the xenon effect in the estimation scheme, whereas shorter windows will imply less xenon variation over time, which discards the xenon dynamics. In [START_REF] Ponçot | Assimilation de données pour la dynamique du xénon dans les coeurs de centrale nucléaire[END_REF], the author use an estimation horizon of two hours with successful xenon estimation.

Input signals applied to the system have faster frequency rates than the xenon dynamics and impact the neutron dynamics. Assuming a discretization of the continuous-time model with frequency T s ∼ 1s would imply N >> 1 to achieve horizons in the range of mere minutes. This is not desirable as the optimization problem would be difficult to solve in a short timespan. Because of the slow dynamics of xenon, such a fast-paced estimator is not necessary anyway.

As such, we have to alter the MHE formulation we presented in the previous section. The estimated states in the estimation window are the result of a numerical simulation instead of a discretized model.

x k+1 = x k + t k+1 t k f (x(τ), u(τ))dτ (4)
y k = h(x k , u k ) . ( 5
)
Since the model is stiff, we use a variable-step numerical integrator between the considered measurement instants and let the integrator choose appropriate time-steps. This implementation imposes that the input signals are to be interpolated at the right time-steps, but ensures a correct input value is used. The integration error is controlled thanks to adaptive schemes, which is essential when long estimation windows are considered.

We assume that the design model is purely deterministic, that is, without process noise. This hypothesis is consistent with the first validation step, which is based on this model. That being said, most of the noise contributions are precisely known on a real installation: the dilution and boration rates are well known, the power extracted from the steam generator is known and the position of the rods is well controlled in normal operation. We also assume that only the beginning of the window is estimated, the remainder of the window can be obtained through integration. Under the previous considerations, MHE is formulated as:

min x[k-N+1|k] J = N i=0 v i 2 R + x[k -N + 1|k] -x k-N+1 2 S k st. x[i + 1|k] = x[i|k] + t k +T s t k f (x(τ), u(τ))dτ v i = y k-i -h( x[k -i|k], u k-i ) x[i] ∈ X, u[i] ∈ U, v[i] ∈ V . (6)
The estimation horizon is a tuning parameter for design. Longer horizons mean that we include more measurements on the system in the cost function and that we are able to capture longer dynamics. The cost function becomes closer to FIE, which is desirable. At the same time, the longer the horizon, the more computation to solve the optimization problem. Horizon length has to be carefully assessed in order to reach satisfactory performances in terms of estimation and time-requirements.

In the current formulation, compared with the "usual" MHE, the designer has an additional parameter to tune: the sampling rate T s = t k+1t k of the measurements included in the cost function. Regardless of the actual measurement rate, we can choose to downsample measurements as to control the complexity of the optimization problem. Horizon length now refers to the quantity (N-1)T s, which corresponds to the actual time spanned by the estimation window, contrary to most implementations for which "N" and "horizon length" are used interchangeably. We will assess the trade-off between higher sampling rate, horizon length and computation time. The load-follow scenario described on figure 2 is used for validation. This scenario is interesting as the successive power variations will disrupt the xenon and iodine dynamics from the equilibrium. Simulation data is representative of actual plant operation as it is obtained in closed-loop using the predictive controller in [START_REF] Dupré | Enhanced flexibility of PWRs (mode a) using an efficient NMPC-based boration/dilution system[END_REF]. We use additive gaussian noise on the measured output, standard deviation for temperature measurements and neutron power being σ t = 0.5 and σ n = 0.2 respectively. Arrival cost is implemented through the constant matrix P. The weigting matrices are :

IV. Results

R = σ t I 2 0 0 σ n I 6 -1 P 1/2 =                    I 12 0 0 0 0 0 10 2 I 9 0 0 0 0 0 10 3 I 2 0 0 0 0 0 10 5 I 6 0 0 0 0 0 I 6                    -1
.

The estimator is initialized with an erroneous estimate x0 with 15% relative error to x 0 , the true initial state. Boron concentration is the only exact initial estimate. Boron absorbs neutrons therefore on the dynamics of the core is similar to xenon. Our tests show that convergence is much slower when starting from an erronated boron concentration. We choose to initialize with a known value, which is a reasonable hypothesis as measurements of the boron concentration are available. Other states could be initialized similarly but do not impact the convergence of the estimator as much as boron.

Relative normalized error (rel) k on the whole state is a weighted norm using P 1/2 so that states have similar contributions to the error, despite their difference in scale [START_REF] Dupre | Design and comparison of two advanced core control sys-tems for flexible operation of pressurized water reactors[END_REF]. For other states we compute the relative error (rel) i,k with respect to the quantity i ∈ {X, I, T, C b , n, c} [START_REF] Lemazurier | Conception d'un système avancé de réacteur PWR flexible par les apports conjoints de l'ingénierie système et de l'automatique[END_REF].

(rel) k = P 1/2 ( x[k|k + N -1] -x k ) P 1/2 x k , (7) (rel)
i,k = xi [k|k + N -1] -x i,k x i,k . (8)
Implementation relied on MATLAB ® [START_REF]MATLAB, Version 9.11.0[END_REF], using fmincon (MATLAB Optimization Toolbox ™ ) with interior-point algorithm to solve the non-linear optimization problem. Numerical integration of the model was handled by the variablestep integrator ode15s. Consumer-grade hardware was used, consisting of an Intel ® Core ™ i5-8365U CPU clocked at 1.6GHz and 16Gb of DDR4 RAM. We provide the average MHE iteration time t iter , which is the average solution time of each optimization problem over the 24h long transient.

A. Trade-off between measurement rate and horizon size

The first validation on the estimation scheme supposes no model mismatch between the design and validation model. We want to assess the tradeoff between horizon length, measurement rate and computation effort. Table I shows the results from running the estimator with different tuning parameters. We see that all estimators are able to reconstruct the state from the noisy measurements (figure 3), with a final relative normalized error (rel) f inal ranging between 0.5% and 1.1%. Focusing on the estimation of xenon, on figure 4, we observe the error decreasing very rapidly in the first iterations of MHE, due to the fact that xenon strongly impacts the dynamics. Experiment 5 yields the best results on relative normalized estimation error. This is to be expected since the horizon in this setting is longer than the others; estimation error at the beginning of the window accumulates in the cost function. When running the same experiments without measurement noise, this difference is even more noticeable. Experiment 4 shows that choosing a "long" horizon is not sufficient, only having half the sampling frequency as experiment 5 while performing significantly worse. We could expect longer horizons to imply significantly longer execution times. It is not the case, showing that the numerical integrator is very effective at choosing the appropriate timesteps, speeding up the evalution of the cost function. The result is a performant MHE execution regardless of the choice of horizon length and sampling period.

B. Sensitivity analysis: a multi-model approach

We assess the robustness of the estimator to model uncertainties. We consider a multi-model approach, in which different models will support MHE, each slighly off from the actual system. As we focus on the xenon/iodine dynamics, we only consider the parameters associated with it. We refer the reader to sections I and II where the parameters are introduced. Table II contains the different model mismatch situations and the error that was introduced on each parameter. Empty entries mean that nominal value was used. The MHE uses the same parameters as experiment 1, as we found that these values produced the best results in this case. (rel) f inal < 2%, while xenon itself is estimated to within an acceptable error level (figure 6): even the worst case scenario (experiment 7) provides X, f inal < 6%. Considering noises on the measurements and the stiff nature of the system, we deem this loss of estimation quality acceptable for real-world application, falling below the noise level of the available measurements.

V. Conclusions

We have designed and implemented a moving horizon estimator for xenon concentration in pressurized water nuclear reactors. The proposed scheme calls for an adapted implementation, taking advantage of variable-step integration. This makes it capable of operating in practise, both because of the quality of the estimate produced, and the computational cost required. Robustness to parameter uncertainties was confirmed through a multi-model sensitivity analysis. This supports suitability of the estimator for realworld applications.

Fig. 2 :

 2 Fig. 2: Reference power throughout the validation transient

Fig. 3 :Fig. 4 :

 34 Fig. 3: Relative normalized state error (in %)

Fig. 5 :

 5 Fig. 5: Relative normalized state error (in %) with model mismatch

Fig. 6 :

 6 Fig. 6: Relative xenon error (in %) with model mismatch

TABLE I :

 I Comparison of different MHE tuning

			Experiment # N	Ts (s) (N -1)T s (s)	t iter (s)
			1	5	300s	1200	7.1
			2	9	150s	1200	8.6
			3	5	150s	800	5.1
			4	5	600s	2400	12.1
			5	9	300s	2400	10.3
		10%					Experiment 1 Experiment 2
	Relative estimation error (%)	1%					Experiment 3 Experiment 4 Experiment 5
		0.1%	0	5	10	15	20
					Time (hours)

TABLE II :

 II Sensitivity analysis of the estimator

			Experiment #	β	K X	Γ X	σ X	t iter (s)
			6				-10%	9.1
			7				-10%	8.4
			8			+5%	6.2
			9			+5%	-10%	6.6
			10		+10%		6.3
		10%					Experiment 6 Experiment 7
	Relative estimation error (%)	1%					Experiment 8 Experiment 9 Experiment 10
		0.1%	0	5		10	15	20
						Time (hours)