Zhuang-Hao Zheng

D.-L Zhang

Bushra Jabar

T.-B Chen
email: chenyx@szu.edu.cn

Mohammad Nisar

Fu Li

Shuo Chen

G.-X Liang

Xiang-Hua Zhang

Ping Fan

Dong-Liang Zhang

Tian-Bao Chen

Yun- Fei Chen

Guang-Xing Liang

Yue-Xing Chen

Realizing high thermoelectric performance in highly (0l0)-textured flexible Cu2Se thin film for wearable energy harvesting

Keywords: Conceptualization, Methodology, Formal analysis, Investigation, Resources, Data Curation, Writing -Original Draft, Project administration, Funding acquisition Thermoelectric, Cu2Se, Flexible thin film, Figure of merit

ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Owing to portable use experience, more and more wearable microelectronic devices have gradually entered the public view [1]. Accordingly, the demand for finding a sustainable self-powered generator for wearable devices instead of using traditionally bulky charging power equipment is urgent [START_REF] Chen | r e -p r o o f Fan[END_REF]. A flexible thermoelectric (TE) thin film device, which can directly converse heat into electricity, is one of the highly expected candidates to realize self-powered by using the temperature difference between the human body and external environment [3]. The efficiency of thermoelectric device is mainly depended on figure of merit ZT (ZT=S 2 σT/κ, where S, σ, κ and T is Seebeck coefficient, electrical conductivity, thermal conductivity, and absolute temperature, respectively.) of the thermoelectric materials. Distinctly, thin films need a high power factor (PF=S 2 σ) and low thermal conductivity to achieve high thermoelectric conversion efficiency in the device. To data, (Bi, Sb)2Te3 based materials are the preferred choice for fabricating thin film device due to their high thermoelectric performance at near room temperature [4,5]. However, compared with equivalent bulks in most reports, the ZT values of (Bi, Sb)2Te3 thin film are low [6,7]. On the other hand, even if high power factor can be produced, large-scale application of (Bi,Sb)2Te3 film will still be limited due to toxicity and high cost of its component; hence, the challenge for thin-film TE device is to find eco-friendly, low-cost, and high ZT value thin films.

In recent decades, eco-friendly bulk thermoelectric materials have made significant progresses as the ZT values had been promoted over 1.5 in some material systems [8][9][10][11][12][13][14]. For example, ZT values reached 2.0 were reported for both PbTe bulk and layered polycrystalline SnSe [9][10]. The GeTe compound via the vacancy suppression and band structure engineering strategy with Bi doing led to a high ZT of 2.0 [11]. However, most of the mentioned materials exhibit the maximum ZT over 600 K, J o u r n a l P r e -p r o o f but have poor TE performance at near room temperature range. Non-toxic and earth-abundant Cu2Se is one of the most studied TE materials which maximum ZT achieved 1.8 for pure Cu2Se [12][13] and the ZT can be further increased to 2.62 by Al doping [14]. The high TE performance is attributed to its special "liquid-like" behavior of the anti-fluorite-structured β-phase. More specially, the experimental result also indicates Cu2Se with monoclinic structured α-phase below 400 K has relatively high ZT value over 0.6 [15] and predicts over 1.1 [16], which is comparable with that of the traditional (Bi, Sb)2Te3. Therefore, Cu2Se provides an exciting commercial opportunity as a potential substitute for fabricating high-efficiency and low-cost flexible thin-film devices. Many works focused on obtaining high-performance Cu2Se thin films have been reported in recent years [17][18][19][20][21]. To data, carrier engineering is the most effective method to improve the TE performance of Cu2Se thin films. For instance, Jin's team proposed that carrier localization increased Seebeck coefficient by fullerenes doping, and the synergistic enhancement of phonon scattering led to reduction of thermal conductivity, resulting in a 20-30% improvement in ZT at 773K [22]. Coincidentally, Chen's group optimized the phonon and carrier transport performance by alloying cuprous selenide with Cu2Se based, achieving a 45% increase in power factor at 300K and an average 24% increase ZT values between 300K and 390K [23]. Carriers transport properties may be further optimized with controllable composition of Cu2Se film, as Cu atoms play an important role in the matrix [24].

In our previous work, a similar strategy was employed which resulted in high power factor of 5.3 μWcm -1 K -2 at 300 K in Cu2Se thin films deposited at solid substrate [24]. The carrier transport of films was significantly improved with enhancement of growth preferred orientation via heat treatment.

Based on the density calculation and the experimental analysis, we confirmed the Cu2Se thin film with (0l0) preferred orientation has higher energy barrier, which will inhibit the bipolar effect and block the J o u r n a l P r e -p r o o f minority carriers, benefiting in achieving high Seebeck coefficient and contributing in the enhancement of power factor. However, Cu2Se thin film prepared at flexible substrate by using the same method has poor electrical transport properties compared with the films deposited on solid substrate due to the poor crystallinity [25]. Therefore, in this work, a facile co-sputtering method was employed to fabricate Cu2Se films on flexible substrate by controlling the accurately composition between Cu and Se, and an annealing process was further used to improve the crystallinity of the films.

Microstructure analysis indicates that thin films have single Cu2Se α-phase with highly (0l0) preferred growth orientation when the Cu/Se composition ratio is over 2.17. These films have high Seebeck coefficient and moderate electrical conductivity, resulting in high power factor. Meanwhile, low thermal conductivity is also found in the Cu-excess thin films due to the introduced point defects in lattice, which can enhance the phonons scattering and lead to low lattice thermal conductivity.

Consequently, a maximum ZT of 0.42 is obtained in the Cu2Se thin films at 275 °C, which is higher than most of flexible Cu2Se based thin films and can even be comparable to the values of bulks at the same measured temperature. Subsequently, traditional silver, molybdenum, aluminium and copper metals are studied as the connecting electrodes of TE Cu2Se legs, and the results show that Ag has the lowest contact resistance. By further introducing a nanometer thickness of Mo buffer layer to prevent the diffusion of Ag, the thermal stability of the as-fabricated device is greatly improved.

Experimental details

A series Cu2Se thin films with different Cu to Se atomic ratios were prepared using the cosputtered Cu target (99.9%) and the Cu2Se alloy target (Cu2Se, ~99.95%). Flexible polyimide (PI, through rate of 70%) with the thickness of 0.125 mm was used as substrate, and cleaned through detergent, ultrapure water, ethanol, and then through the compression air jet equipment and drying box J o u r n a l P r e -p r o o f for drying treatment. The target powers of Cu2Se target and Cu target were selected to 55 W and 15 W, respectively. The background pressure was 8×10 -4 Pa, and working pressure was adjusted to 0.5 Pa by setting the argon gas of 40 sccm. The deposition time of Cu2Se is 30 min. The samples were annealed in a glove-box heating system in a nitrogen atmosphere at 300 ℃ for 30 min. The thickness of the thin films and composition measured by scanning electron microscope (SEM) and energy dispersive X-ray spectroscopy system (EDS). The thickness of the film was measured using a profilometer (Dektak XT, Bruker). In this work, the atomic ratios of Cu/Se of films are increase gradually with increasing the operating time from 0 min to 11 min ( 0, 3, 5, 7, 9 and 11 min, respectively) of Cu target during the co-sputtering process. The corresponding samples are named as Cu-0min, Cu-3mins, Cu-5mins, Cu-7mins, Cu-9mins and Cu-11mins, respectively. The corresponding atomic ratios of Cu/Se of films are 1.99, 2.05, 2.17, 2.20, 2.26 and 2.31, respectively. The EDS results are listed in Table S1 (supporting information).

The crystal structure of the film was analyzed by X-ray diffractometer (XRD, Ultima-iv) and the test Angle ranged from 10°to 60°with the scan speed is 2°/min with the CuKα radiation. Preferred orientation information of Cu2Se thin films was characterized by four-circle X-ray diffractometer (X'PERT, Philips) using pole figure scans. The surface morphology was measured by scanning electron microscope (SEM,Supra55), and depth microstructure was also measured by transmission electron microscopy (TEM, JEM3200FS). The chemical valence states of copper and selenium were characterized by X-ray photoelectron spectroscopy (XPS, Escalab 250Xi) using a monochrome Al-Ka X-ray source. Plasma etching was performed for 30 seconds to remove surface impurities. The electronic band structure and density of states (DOS) for Cu2Se and Cu-excess Cu2Se are performed by first-principles plane-wave pseudopotential formalism as implemented in the VASP package, the J o u r n a l P r e -p r o o f calculation detail can be found in Supporting Information.The thickness was measured by a step meter and scanning electron microscope. Carrier concentration was provided by the Hall tester (HL5500PC, Nanometrics) under a reversible magnetic field (0.5 T) by the van der Pauw method. Electrical conductivity and Seebeck coefficients as a function of temperature were provided by the Seebeck conductivity coefficient measurement system (SBA458, Netzsch). Thermal conductivity was measured using the Van-der-Pauw test [26]. The method strictly requires the film to be covered on the chip contact, and the transverse thermal conductivity is measured after the probe touches the contact.

The film preparation can be achieved using a pull-0ff foil deposition mask, which is suitable for most film-making processes. The measurement was performed from room temperature to 275 ℃. Keithley 2400 was used to measure IV data of thin film devices and plot curves. The combined uncertainty for the experimental determination of ZT is about 15% and is caused by the three respective measurement parameters including σ, S and κ.

Results and discussion

To detect structural changes and understand the variation in lattice of Cu2Se with the change of stoichiometric ratio, we naturally carried out X-ray diffraction analysis. Fig. 1(a) shows XRD patterns of Cu2Se based films with different stoichiometric ratios at fixed heat treatment temperature. The films presents high temperature phase (β-Cu2Se, F-43m) for Cu1.99Se and Cu2.05Se samples, which gradually convert to low temperature phase (α-Cu2Se, Fm-3m) as the increase of Cu content. It is worth noting that high temperature phase β-Cu2Se is stable at room temperature, which may arise due to the existence of low dimensional effect in the thin film, consistent with the results of the previous study [18]. For the films with composition of Cu-7mins, Cu-9mins, and Cu-11min, the XRD patterns are well indexed as α-Cu2Se phase (PDF#37-1187), and no second phase is observed. Furthermore, the Rotating counterclockwise around the vertical X-axis is 0° to 360° in polar coordinates. It can be seen that the test results at all locations indicate that the film has peak intensity gathering points at Phi=35.71°, Phi=95.24° and Phi=152.53°, which exactly represent the three planes of (030), ( 060) and (090), consisting with the (0l0) preferred orientation in XRD patterns (Fig. 1(a)) As shown in Fig. S3 (Supporting Information), there are other peak densities in the film, but their intensity is not as high as the above three positions, so it does not affect the judgment of preference orientation. XPS measurements were employed on Cu2Se based thin films (Cu-3mins, Cu-5mins, and Cu-7mins) to study the binding states. As shown in Fig. 1(g), the core-level spectrum exhibits two strongest peaks at ~932.4 eV and ~952.2 eV, respectively, namely, Cu 2p3/2 and Cu 2p1/2. It basically matches the spin-orbit phenomenon of Cu (I) [27]. More notably, we do not observe impurity peaks in the spectra, excluding the existence of copper oxidation in this sample. The characterization of Se (-Ⅱ) in a consistent, cohesive environment [28] is described intuitively in Fig 1 (h). It can be seen that the nearly symmetrical wide peak spans from 54 to 56 eV and this wide peak can be divided into two small peaks, which represent Se 3d5/2 and Se 3d3/2, respectively, and refer to two positions ~54.2 and ~54.9 eV. To sum up, the chemical states of elements in thin films prove that only Cu 1+ and Se 2-bonding states can be observed in as-fabricated films. In addition, the weak dependence of peak positions for The surface of morphology of the Cu2Se-based thin film is studied by SEM, as shown in Fig. 2.

It can be found that Cu-0min (Fig. 2a), Cu-3mins (Fig. 2c), as well as Cu-5mins (Fig. S4a) thin films J o u r n a l P r e -p r o o f are growing dense which grains are uniformly distributed with the size ~ 50 nm. The cross-section images reveal the thickness of films is ~400 nm, and the grains are irregular arranged (Fig. 2b). As for the films with higher Cu content, namely Cu-7mins (Fig. 2d), Cu-9mins (Fig. 2e) and Cu-11mins (Fig. S4b, Supporting Information), most of the grains are growing vertical to the surfaces, which can be further confirmed via the cross-section images (Fig. 2e). The regular growth of the grains has consisted with the high orientation feature as concluded from the XRD patterns and polar diagram analysis (Fig. 1). Fig. 2 (g ~ i) show element mappings of Cu-7mins film, both Cu and Se elements are homogeneously distributed, and no regional aggregation can be observed according to the detectable limit of EDS. Hence, the existence of the extra Cu in the films will be further explored via the microstructure analysis. [1] in (b), the insets is IFFT images of white-dotted region in (c)), (d) enlarged HRTEM image of selected area [START_REF] Chen | r e -p r o o f Fan[END_REF] in (b), (e) HRTEM image of grain boundary between two crystalline phases clearly indicates coherent phase boundary (selected area [3] in (b)), (f) IFFT images of an area from a crystalline phase.

The extremely high (0l0) orientation and the significant improvement of TE performance are J o u r n a l P r e -p r o o f closely related to the excess Cu content for Cu2Se film; hence, a transmission electron microscopy (TEM) was performed to investigate the microstructure of Cu-7mins film to confirm the existence of excess Cu. Fig. 3 (a) shows the low magnification TEM image of Cu-7mins film, in which the highresolution TEM image marked with white-dotted rectangle is presented in Fig. 3 (b). It can be seen that the grains are embedded in the Cu-7mins film. Fig. 3 (c) and Fig. 3 (d) show enlarged HRTEM images of selected regions [1] and [START_REF] Chen | r e -p r o o f Fan[END_REF] in Fig ( b). The images clearly show the semi-coherent boundary between the highly crystalline phase and the polycrystalline phase. Of course, the crystal plane spacing of (012) in the two figures is approximately 2.062 Å, agreeing with the matrix of Cu2Se. The HRTEM image (Fig. 3(e)) of grain boundary between two crystalline phases indicates coherent phase boundary.

The inverse fast Fourier transform (IFFT) image of an area from a crystalline phase is shown in Fig. S1, Supporting Information). Subsequently, the presence of highly crystalline phase and the polycrystalline phase illustrate the atomic models of the two orientations, agreeing with the high (0l0) orientation of as fabricated film. In addition, the stress arising from the semi-coherent interfaces in the film (Fig. S5 (Supporting Information)) is conducive to enhancement of phonon scattering and leads to the reduction of the thermal conductivity. [19]. To explain the variation of and S, carrier concentration n was measured at room temperature, and the carrier mobility μ was calculated J o u r n a l P r e -p r o o f from the experimental  and n ( = neμ, where e is the carrier charge) With increasing Cu content from 1.99 to 2.31 in Cu2Se films, n is gradually reduced, which is the major factor for the enhancement of S and decrease of . That is because the excess Cu will reduce the Cu vacancies and results in the decline of n, consisting with the experimental results of CuxSe bulks reported by Shi et al [29].

3(f). The IFFT image represents the atomic distribution within α-Cu2Se atomic framework (Table

Meanwhile, the mobility of α-Cu2Se films increased with Cu content (Cu/Se= 2.17~2.31) (Fig. 4 (c)), corresponding with the reduction in calculated effect mass m* (Fig. S6a-6b, Supporting Information), which was calculated using the single parabola band (SPB) model. Obviously, relative high m* of films with excess Cu is responsible for the decent S. In addition, due to the phenomenon of high preferred orientation of α-Cu2Se films, the carrier concentration and mobility of the film are generally higher than that of other reports [20,30]. The band structures of β-Cu2Se (Fig. 4(d)) and α-Cu2Se (Fig. 4(e)), and corresponding density of states (DOS) (Fig. S7(a)-(b), Supporting information) exhibit narrow band-gap semiconductor characteristic with different orbital energy distribution near the valence band maximum (VBM), which is believed to be caused by the disorder of Cu atoms [31]. It is noteworthy that the slope of density of state (DOS) of α-Cu2Se near VBM is steeper than that of β-Cu2Se, consisting with the experimental results of higher Seebeck coefficient of α-Cu2Se (Fig. 4(b)).

Due to the regularly variation in  and S, a maximum power factor value ~ 8.4 μWcm -1 K -2 at room temperature for Cu2.20Se thin film, and then increase to 11.12 μWcm -1 K -2 at 275℃ (Fig. 4 (f)). As shown in Fig. 4 (g), the high PF is higher than that of other reported films deposited by different methods [17-20, 25, 32-34], including sputtering [17,32], thermal evaporation [18], pulsed laser deposition [19], solution immersion method [20], sputtering combined with thermal evaporation method [25], ion beam sputtering [33], spin coating process [34]. The high PF is also comparable with that of Cu2Se bulk (PF=12.52 μWcm -1 K -2 ) reported by Yu et al.. [35] With regard to the thermal J o u r n a l P r e -p r o o f transport properties of the as-fabricated Cu2Se films, the thermal conductivity  was measured along the in-plane direction (same measured direction with electrical transport properties), the temperature dependence of  is plotted in Fig. 4 (h). With increased content of Cu, the  and l (l is determined from different value of  and electrical, thermal conductivity e, seen the Fig. S8, Supporting Information) are significantly reduced. This may be attributed to the significantly enhancement of phonons scattering in the films, as the excess Cu atoms can be easily formed as phonon scattering centers and the existence of stress arising from the semi-coherent interfaces helps to block the transport of phonons. Finally, the ZTs of the Cu2Se films were determined (Fig. 4 (i)). A ZT value of 0.21 is achieved for Cu-7mins film with Cu:Se atomic ratios of 2.20 : 1 at room temperature, which value gradually increase to 0.42 at 275 ℃. Based on the good thermoelectric properties as demonstrated above, we further investigate the practical applicability of the Cu2Se based flexible thin films. It is well-known that the total resistance of the device should be as low as possible to ensure high performance of the device. Therefore, selecting a suitable and accurate metal electrode is one of the key procedures to assembly device as the resistance of composed material itself is almost unchanged during the operation. We, firstly, deposited different kinds of metal electrodes on as-fabricated Cu2Se films and measure U-I characteristic curve to check the contact resistance between the film and metal electrodes (Fig. 5a).

J o u r n a l P r e -p r o o f

The contact between metal and Cu2Se film is divided into Schottky contact and low resistance ohmic contact, the actual situation of which kind of contact depends on the material characteristics between metal and semi-conductor [36][37][38]. Four common single metal electrodes, including Cu, Al, Mo and Ag are performed in this work. Fig. 5 (a) to Fig. 5 (c) show the basis tests for electrode selection, based on TLM (Transmission Line Model) measurement principle [39], which model construction and physical demonstration is provided in Fig. S9 (Supporting Information). TLM is a measurement method that uses metal electrodes of the same length and width to be uniformly arranged on the thin film. The total resistance between adjacent metal electrodes can be measured by the difference of electrode spacing. The IV curves of each metal electrode (Fig. 5(a)) show the highest voltage when J o u r n a l P r e -p r o o f employing Ag electrode under the same current. This phenomenon also precisely corresponds to the lowest total resistance (RT) obtained by different metal electrodes at different distances, as shown in Fig. 5 (b). It is worth mentioning that the theoretical value of resistance between two adjacent electrodes can be calculated by the formula [40]:𝑅 𝑇 = 2𝑅 𝑚 + 2𝑅 𝑐 + 𝑅 𝑠𝑒𝑚𝑖 , where Rm is the resistance of the metal electrode, which is usually very small and can be ignored; RC is the contact resistance between the metal and the semiconductor, Rsemi is the resistance of the semiconductor. The actual situation is that RC often changes with temperature, that is, the temperature difference will bring about the difference between theory and practice. In order to facilitate calculation, RC is a fixed value by default.

Therefore, the difference between theoretical resistance and experimental resistance will not be large and within a controllable range. We also check the contact resistivity (ρ c ) to characterize the ohmic contact quality, the ρ c is generally expressed as following formula: 𝜌 𝑐 = 𝑅 𝐶 𝐿 𝑇 𝑊 where LT is the spacing between electrodes and W is the width of single electrode. The lower the contact resistivity means the better ohmic contact quality. It can be found that Ag electrode exhibit lower ρ c (0.36×10 -2 Ω•cm 2 ) than other metal electrodes. However, Ag electrode is not adopted as the final choice as the Ag diffuses into Cu2Se layer (Fig. 5d) after heat treatment above 275℃ (the cross-sections of Ag electrode bonded with Cu2Se before heating can be found in Fig. S10a,Supporting Information), which deteriorates the chemical stability of as fabricated film. The deterioration in ohmic contact can be further confirmed from the result of U-I curve, RT and ρ c in Fig. 5a-5c. Hence, we set a Mo layer with a thickness of ~ 50 nm, which possesses better chemical stability, between Ag electrode and Cu2Se based film (See the SEM image in Fig. 5(e)), expecting to block the diffusion of Ag. Qiu et al. also used the Mo as diffusion barriers on the Cu2Se/Yb0.3Co4Sb12 TE couples used in the module to block the the diffusion of Cu [41]. The experimental results also proves that this approach is effective. The J o u r n a l P r e -p r o o f films bonded with composited Ag/Mo electrode shows a much higher voltage than the other metal electrodes at the same current (Fig. 5 (a)), as well as a lower total resistance (Fig. 5 (b)) and better contact quality (Fig. 5 (c)). More importantly, the Ag diffusion is effectively restrained after heat treatment by using the composited electrode, which is clearly shown in Fig. 5 (e). The cross-section mapping images of Ag/Mo composited electrode bonded Cu2Se film before heating can be found in

Conclusions

In this work, highly (0l0) orientation flexible Cu2Se films with variable Cu/Se ratios were fabricated by a facile co-sputtering method. Excess Cu in the Cu2Se films played important roles in optimizing the carrier concentration by suppressing the formation of Cu vacancies and donating electrons, leading to the high Seebeck coefficient. The highly (0l0) preferred orientation contributed to high carrier mobility and maintained decent electrical conductivity. As a result, an excellent PF of 11.12 μWcm -1 K -2 was obtained for Cu2Se based film at 275℃. Combined with the low thermal conductivity of as-fabricated films, a maximum ZT of 0.42 was achieved at 275℃. Subsequently, the thermal stability of the flexible thermoelectric device assembled with high performance Cu2Se films was greatly improved after introducing a Mo buffer layer into electrode layer. The maximum power density of 4.28 Wm -2 was achieved at a temperature difference of 50 ℃, demonstrating its potential in wearable applications.

  Fig. 1(c) and Fig. 1(d). It can be analyzed that there are high peak intensities at 2θ=13.039° (030) and

  J o u r n a l P r e -p r o o fCu and Se elements indicates the stable valence states for all films.

Fig. 1 .

 1 Fig.1. (a) XRD patterns of thin films with different composition ratios at the same heat treatment

Fig. 2 .

 2 Fig.2. Surface morphologies of the film with different Cu:Se atomic ratios of (a) 1.99:1 (Cu-0min), (c)

Fig. 3 .

 3 Fig.3. TEM analysis of Cu-7mins thin film; (a) low resolution TEM image, (b) high resolution TEM

Fig. 4 (

 4 Fig.4(a) and Fig.4(b) show the temperature dependence of electrical conductivity  and Seebeck

Fig. 4 .

 4 Fig.4. (a) Electrical conductivity of Cu2Se-based films at different stoichiometric ratios. (b) Seebeck

Fig. 5 .

 5 Fig.5. Comparison of (a) IV curves, (b) total resistance RT and (c) contact resistivity ρ c when using

Fig. S10b (

 S10b Fig.S10b (Supporting Information). A flexible thermoelectric device consisting of 10 Cu2Se films (26

Acknowledgements

This work is supported by National Natural Science Foundation of China (No.11604212), Guangdong Basic and Applied Basic Research Foundation (2019A1515110107 and 2020A1515010515), and Science and Technology plan project of Shenzhen (20200811230408001).

The authors are thankful for the assistance on HAADF-STEM observation received from the Electron Microscope Center of the Shenzhen University.

Highlights

(1) Flexible Cu2Se based thin films with desired highly (0l0)-preferred orientation were successfully achieved by optimized growth temperature and composition;

(2) Introduced excess Cu in flexible Cu2Se films contribute to the improvement electrical and thermal transport properties simultaneously;

(3) A Mo buffer layer is innovatively introduced into electrode layer, greatly enhance the thermal stability of flexible thermoelectric device assembled with high performance Cu2Se films.

J o u r n a l P r e -p r o o f

Declaration of interests

☒ The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

☐ The authors declare the following financial interests/personal relationships which may be considered as potential competing interests:

J o u r n a l P r e -p r o o f


Fabien Robin

MD MSc Ludivine Dion

MD PhD Vincent Lavoue

MD MSc Anna Goujon

MD Karim Boudjema

MD Laurent Sulpice

MD Xavier Morandi

Dr Nyangoh Krystel

Timoh

Robin Fabien

MD Krystel Nyangoh Timoh
email: krystel.nyangoh.timoh@chu-rennes.fr

Relationship between uterine veins, ureter, and hypogastric nerve for uterine transplantation: An anatomic study

Keywords: Uterine transplantation, anatomy, uterine veins, ureter, hypogastric plexus

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

INTRODUCTION

Uterus transplantation (UT) is the only available treatment for women with absolute uterine factor infertility (AUFI) caused by congenital or acquired absence of the uterus. In this setting, Brännström's team demonstrated the feasibility of the procedure for women who desired to become pregnant, and announced the first birth in October 2014 [START_REF] Brännström | First clinical uterus transplantation trial: a six-month report[END_REF][START_REF] Brännström | Livebirth after uterus transplantation[END_REF]. To date, more than 76 transplants have been performed worldwide, resulting in the births of 19 live children. UT is therefore a growing and promising technique for women with AUFI to achieve motherhood [START_REF] Dion | La transplantation utérine : essor de la technique chirurgicale, nouvelles indications. État de l'art avec revue de la littérature[END_REF]. As most grafts are obtained from living donors (> 90% of UT performed), the main issue in the procedure is to prevent complications for the living donor while respecting the integrity of the graft. Two veins can be selected for removal: the utero-ovarian vein, which requires oophorectomy in the donor, and the uterine vein. According to Brännström, the most challenging step is the removal of the uterine veins, which is associated with an elevated risk of: ureteral injury (estimated at over 20%) [START_REF] Chmel | Robot-assisted donor hysterectomy in uterus transplantation -a modality to increase reproducibility[END_REF]; venous lesions and subsequent thrombosis; and a greater operating time [START_REF] Brännström | First clinical uterus transplantation trial: a six-month report[END_REF]. Furthermore, cases of postoperative bladder dysfunction in the donor could be related to this operating step [START_REF] Chmel | Revaluation and lessons learned from the first 9 cases of a Czech uterus transplantation trial: Four deceased donor and 5 living donor uterus transplantations[END_REF]. However, in contrast to knowledge of the anatomy of the hepatic pedicle in the field of liver transplant [START_REF] Hiatt | Surgical anatomy of the hepatic arteries in 1000 cases[END_REF], little is known about the anatomy of the uterine pedicle and the uterine veins. In classic textbooks, two uterine veins are described within the uterine pedicle [START_REF] Gray | Grays Anatomy[END_REF], but the intimacy of their relationships to the ureter and the uterine artery within the parametrium is unclear. Also, there is a lack of knowledge of their caliber and of where they terminate in the internal iliac vein. These facts are as crucial for deceased-donor uterine procurement as for live-donor retrieval. Finally, available descriptions of uterine vein removal during uterine procurement are succinct.

We therefore aimed to provide a comprehensive standardized description of the venous drainage of the uterus in relation to the ureter and pelvis nervous tract through investigation of deceased subjects and digital reconstruction of female pelves, with a view to establishing surgical anatomical landmarks for both living donor procurement and transplantation.

METHODS

Study design

The anatomical data were obtained from two types of source: i. sixteen hemipelves from eight female anatomical subjects dissected at the Rennes School of Surgery (anatomy laboratory of Rennes University); ii.digital reconstruction with the Anatomage® Table of four hemipelves from two additional female subjects.

The scientific committee of Rennes University ensured that written consent for body donation had been obtained and filed prior to death for all the anatomical subjects.

Our work complied with French regulations. The study was exempt from French law pertaining to biomedical research (the Huriet-Serusclat Law, 20 December 1998) as no additional interventions were required.

Anatomical study

Terminology

The second edition of the Terminologia Anatomica was used to avoid confusion in describing anatomical terms for the pelvic connective tissues [START_REF] Fipat | Terminologia Anatomica[END_REF]. When a term was missing in the Terminologia Anatomica, we used unofficial but frequently cited terms [START_REF] Ercoli | Terminologia Anatomica versus unofficial descriptions and nomenclature of the fasciae and ligaments of the female pelvis: A dissection-based comparative study[END_REF].

Procedure

The dissection protocol consisted of making a median longitudinal abdominal incision from the xiphoid to the median of the pubic symphysis. For better exposure of the pelvis, the sigmoid colon and the upper rectum were sectioned and the entire digestive tract was reclined upwards. In four of the sixteen subjects, the dissection was done after arterial casting with red colored latex.

The following steps were repeated during the dissection of the 16 female hemipelves: 1) Identification of the retro-ligamentous part of the ureter.

We made a peritoneal incision to open the pararectal fossa at the lateral wall of the pelvis at the level of the iliac artery bifurcation.

2) Identification of the internal iliac vein.

We opened the medial pararectal space, also called the Okabayashi space, backward and laterally to the uterosacral ligament. The hypogastric nerve was identified within this space.

We then identified a vein that passed below the ureter, which we called the deep uterine vein (DUV), within the paracervix.

We also identified the inferior hypogastric plexus (IHP) and the endings of the pelvic splanchnic nerves (PSN) within the paracervix.

3) Identification of the ligamentous part of the ureter.

We dissected the ureter carefully downwards until it tunneled through the parametrium and crossed the ureterine artery, determining the ligamentous part of the ureter. We respected the meso-ureter, the anatomical structure that corresponds to a thin, sagittally-oriented connective mesentery enveloping the ureter within the presacral fascia [START_REF] Ercoli | Terminologia Anatomica versus unofficial descriptions and nomenclature of the fasciae and ligaments of the female pelvis: A dissection-based comparative study[END_REF].

We then identified the ligamentous part of the uterine artery. 4) Identification of a vein that passes above the ureter, which we called the superficial uterine vein (SUV), within the parametrium.

5) When possible, the vaginal artery was identified.

6) The uterine veins were followed until their termination in the internal iliac vein.

7) Identification of the venous plexus.

We sectioned the round ligament and then opened the anterior leaf of the broad ligament to access the mesometrium. Here, we identified the venous plexus and the ascending branch of the uterine artery.

Digital reconstruction with the Anatomage Table

The Anatomage Table (Anatomage, Inc., San Jose, CA) is a virtual dissection table created through imaging of real human anatomical subjects. Reconstructions combine data from the Visible Korean Project and the Visible Human Project. In both projects, cadavers were subjected to MRI and CT scanning and were then cut into sections ranging from 0.2 to 1 mm using a cryomicrotome; the slices were photographed with a digital camera. By combining computerized data, the table produces a real cadaver bed and demonstrates real patient data on a life-size scale [START_REF] Ackerman | The Visible Human Project: a resource for education[END_REF][START_REF] Park | Visible Korean Human: its techniques and applications[END_REF][START_REF] Martín | Possibilities for the use of Anatomage (the Anatomical Real Body-Size Table) for Teaching and Learning Anatomy with the Students[END_REF]. These data are used for research with the agreement of the company Anatomage Inc. We provide a digital reconstruction of the pelvic area from a deceased subject to illustrate our data highlighting the uterine veins, uterine artery, ureter, and the inferior hypogastric plexus.

Results

Results from dissection of anatomical subjects.

Uterine veins

In all 16 hemipelves we identified the uterine venous plexus in the broad ligament within the mesometrium. In all cases, the venous plexus was accompanied by the ascending branch of the uterine artery (Figure 1). From this uterine venous plexus, one or two uterine veins emerged.

The superficial uterine vein

In four of the 16 hemipelves we identified the SUV coming from the uterine plexus and passing above the ureter. The SUV was very thin and always accompanied by the uterine artery at its caudal aspect. Thus, the SUV is within the parametrium (Figures 2 and3).

The deep uterine vein

In all 16 hemipelves, a DUV was identified coming from the uterine plexus and passing below the ureter. The DUV has a caliber double or triple that of the SUV and is located within the paracervix. It is not a satellite to the uterine artery and runs posteriorly, cranially and laterally until it joins a short segment in common with the SUV. The two veins finally anastomose into a common uterine vein (Figure 3), which then terminates in the internal iliac vein at various levels.

Relationship with the pelvic ureter

The pelvic ureter closely related to the uterine veins.

Within the parametrium, the ligamentary ureter is crossed by the SUV (when present) and the uterine artery, and the DUV passes underneath. Thus, the ligamentary ureter is surrounded by uterine veins -the SUV in a cranial position, and the DUV in a caudal position -forming a triangle into which the ureter descends. The base of the triangle corresponds to the lateral border of the uterus, the upper side the SUV, and the lower side the DUV. The whole structure is directed laterally (Figures 2 and4).

Relationship of the uterine veins to the pelvic autonomic nervous system

In all 16 hemipelves, the inferior hypogastric plexus (IHP) and pelvic splanchnic nerves (PSN) were identified in the paracervix. They were always situated below the ureter and the DUV. Thus, in the lateral aspect of the uterus, cranial-to-caudal, lie the uterine artery, the SUV (when present), the ureter, and the DUV. Consequently, the DUV is the most cranial [caudal?] element of the paracervix forming the roof of the IHP and PSN. Finally, the DUV is positioned between the ureter cranially and the IPH caudally.

Results from virtual dissection

The digital reconstruction of the pelvic area illustrates the results (Figure 5). The SUV was not seen in either of the two subjects. The DUV was located below the ureter and the IHP and PSN below the DUV. We made an animated movie of these results (Annex 1).

.

DISCUSSION

The present study provides a comprehensive description of the anatomy of the uterine veins, and highlights crucial relationships. We report that: (1) there are two uterine veins -the small and inconstant SUV, which runs cranially to the ureter, and the larger DUV, which runs caudally to the ureter; (2) these two veins describe an anatomical triangle with the medial uterine border where the ureter passes; and (3) the DUV is "sandwiched" between the ureter cranially and the hypogastric nerve plexus caudally. These findings have major surgical and clinical relevance to the field of UT.

In gynecology, while the uterine artery is precisely described, the anatomy of the uterine veins remains unclear. Moreover, definitions and description of the uterine veins are missing in the Terminologia Anatomica [START_REF] Fipat | Terminologia Anatomica[END_REF]. The most detailed description of the pelvic veins seems to have been provided in 2004 by Kamina, who highlighted their avalvular character and their topography in the parametrium and paracervix [START_REF] Kamina | Anatomie Clinique[END_REF]. However, unlike the arteries, which have been extensively studied and are assessed before transplantation [START_REF] Johannesson | Vascular pedicle lengths after hysterectomy: toward future human uterus transplantation[END_REF][START_REF] Orhan | From Where Does the Uterine Artery Originate? A Prospective, Observational Laparoscopic Anatomic Study[END_REF][START_REF] Leonhardt | Imaging evaluation of uterine arteries in potential living donors for uterus transplantation: a comparative study of MRA, CTA, and DSA[END_REF], the uterine vein dissection step is not well described. This is in contrast to liver transplantation, where the different types and frequencies of vascular anatomical variations are now well known, allowing the surgeon to anticipate difficulties and determine a reimplantation strategy [START_REF] Hiatt | Surgical anatomy of the hepatic arteries in 1000 cases[END_REF]. Additionally, to the best of our knowledge, there is a lack of original anatomical work on the anatomy of the uterine veins for UT. However, the step of uterine vein dissection is associated with pre-and post-operative complications for both the living donor and the graft [START_REF] Jones | Human uterine transplantation: a review of outcomes from the first 45 cases[END_REF]. The dual digital and cadaveric approach developed in this work is innovative. Anatomage® data are derived from a combination of high-definition imaging data (CT, MRI) and cadaveric sections and thus constitute real high-fidelity anatomical models rather than a theoretical vision. This allowed us: (1) to evaluate the objective anatomical relationships of the veins in situ without the structural displacement secondary to surgical dissection artefacts; (2) to confirm the results observed in dissection, thus reinforcing the main message of the study; and (3) to illustrate the results for a better understanding of the uterine veins.

In the most recent series of uterine living donors, the rate of postoperative grade IIIB complications was over 11% [START_REF] Dindo | Classification of Surgical Complications[END_REF][START_REF] Jones | Human uterine transplantation: a review of outcomes from the first 45 cases[END_REF]. This complication rate is similar to those reported for liver (10%) and kidney (8%) living donors. Among the major complications, ureteral injuries are predominant (20%) and far more frequent than in classic intra-or extra-fascial hysterectomies [START_REF] Timoh | Standardization of laparoscopic extrafascial hysterectomy: anatomic considerations to protect the ureter[END_REF][START_REF] Testa | The Evolution of Transplantation From Saving Lives to Fertility Treatment: DUETS (Dallas UtErus Transplant Study)[END_REF].

Ureteral injuries negatively affect quality of life, potentially requiring a secondary ureter reimplantation in the living donor [START_REF] Fageeh | Transplantation of the human uterus[END_REF][START_REF] Brännström | First clinical uterus transplantation trial: a six-month report[END_REF]. In the present study, we emphasize that the ureter is located in the center of a utero-venous triangle formed by the SUV cranially, the DUV caudally, and the lateral aspect of the uterus medially.

Knowledge of this anatomical landmark is critical for protecting the ureter.

The success of other solid organ transplantations, such as the liver, has been achieved through precise knowledge of the relevant vascular anatomy and its variations. While the problem of inflow through the uterine artery is not an issue, that of venous outflow is a subject of debate, and routine practice varies across the world. Double ovarian and uterine drainage theoretically allows two options for venous anastomosis. Exclusive ovarian venous drainage would seem interesting because it avoids the need to dissect the uterine veins within the parametrium. However, it requires oophorectomy in 78% of cases [START_REF] Dion | La transplantation utérine : essor de la technique chirurgicale, nouvelles indications. État de l'art avec revue de la littérature[END_REF] with the risk of induced menopause, and also has the highest failure rate of UT procedures (38%) [START_REF] Jones | Human uterine transplantation: a review of outcomes from the first 45 cases[END_REF]. In addition, the ovarian vessels are of good caliber only in young donors.

In the coming years, with an increase in older donors, these veins could be too small.

Thrombosis is a threat in the event of poor drainage but can be avoided by good caliber regular veins and integrity of the vessels. We describe a duality of uterine veins -the DUV and the SUV -which until now has been poorly identified. Of the two, the DUV represented the larger and more constant drainage vessel in our dissections and we would therefore recommend this vein for graft drainage. The SUV, which was observed in 25% of our dissections, has probably been underestimated because of its fragility and small size in anatomical subjects. If both uterine veins are present, one solution could be to cut the SUV, procure the graft without injuring the ureter, and re-anastomose the transected SUV at the back table [START_REF] Brännström | Livebirth after uterus transplantation[END_REF].

In their review of the first 45 procedures, Jones et al. reported 28% of grade I-II complications, including consequences of nerve injury [START_REF] Jones | Human uterine transplantation: a review of outcomes from the first 45 cases[END_REF]. These included postoperative functional sequelae such as fecal incontinence and urinary post-voiding dysfunction. The Czech team accurately described how to resolve bladder hypotonia by prolonged suprapubic catheterization [START_REF] Chmel | Revaluation and lessons learned from the first 9 cases of a Czech uterus transplantation trial: Four deceased donor and 5 living donor uterus transplantations[END_REF]. Our work provides an explanation for these complications by describing the close caudal relationship of the DUV to the inferior hypogastric plexus. A nerve-sparing technique was standardized for endometriosis surgery by the team of Ceccaroni, who described the DUV as an anatomical landmark for nerve-sparing [START_REF] Ceccaroni | Laparoscopic nerve-sparing transperitoneal approach for endometriosis infiltrating the pelvic wall and somatic nerves: anatomical considerations and surgical technique[END_REF]. This technique decreases the risk of postoperative voiding dysfunction. However, even if the autonomic system can be identified, reproducible, nervesparing techniques have not been developed for all gynecological pelvic surgeries. We describe here the "sandwich" situation of the DUV with the ureter above and the IHP below as a major landmark for uterus harvesting.

On the basis of these results, with a view to protecting both the ureter and the autonomic system during uterine vein dissection, we propose the following: ( 1) to open the pararectal space to identify the IHP and PSN for nerve-sparing, (2) to identify the DUV, which is clearly visible within the paracervix owing to its greater caliber, (3) to dissect the retro-ligamentary ureter until it crosses the uterine artery, and (4) to identify the SUV. Here, the choice of sacrificing the SUV is possible. Once the DUV is free from the ureter cranially and the nerves caudally it can be removed safely.

The development of robotic pelvic surgery in recent years has led to increased accuracy in dissection. In both male and female surgery, nerve-sparing has become easier with improved vision [START_REF] Magrina | Robotic nerve-sparing radical hysterectomy: feasibility and technique[END_REF]. This approach is being developed to minimize complications and trauma in the living donor and a trial evaluating the procedure is ongoing (NCT02987023) [START_REF] Brännström | Evolution of surgical steps in robotics-assisted donor surgery for uterus transplantation: results of the eight cases in the Swedish trial[END_REF][START_REF] Johannesson | Dallas UtErus Transplant Study: Early Outcomes and Complications of Robot-assisted Hysterectomy for Living Uterus Donors[END_REF]. In addition, endometriosis surgery has emphasized the need to improve knowledge of the extraperitoneal region, particularly that of the nerves [START_REF] Fermaut | Deep infiltrating endometriosis surgical management and pelvic nerves injury[END_REF][START_REF] Kanno | Robot-assisted Nerve Plane-sparing Eradication of Deep Endometriosis with Double-bipolar Method[END_REF]. Augmented reality techniques and the use of indocyanine green should increase the surgeon's ability to discriminate structures.

Despite our results, this study has some limitations. First, we included only a small number of anatomical subjects and, consequently, we probably did not capture the full range of variations. For example, termination of the DUV in the pudendal or gluteal veins can be observed in gynecological surgery, although it has not been described in the literature. Future studies on living donors are required to focus more on the modal and non-modal uterine vein anatomy to establish the full vascular pattern. Second, we worked mainly on elderly and menopausal cadaveric donors. Advanced age is accompanied by involution of the internal genital organs leading to a moderately different anatomy from that of the female candidate for uterine procurement. Finally, the small size and vacuity of the SUV could have led to their section during dissection, thus underestimating their frequency.

Conclusion

This study establishes a standardized anatomy of the uterine veins. It describes the very close relationship of the ureter passing between the uterus and SUV in the parametrium but also the posteroinferior relationship with the IHP. These data provide additional knowledge to the surgeon performing UT both to limit donor injury and to choose the best graft drainage. In green: the urinary system: the bladder, urethra and the ureter. In purple: the uterus and the vagina. In yellow, the inferior hypogastric plexus.

FiguresFigure 1 :

 1 Figures 2B and 4 were reproduced from « Rubod C, Collinet P. Chirurgie de l'endométriose -Approche multidisciplinaire médico-chirurgicale. Paris 2022 Elsevier Masson SAS. Tous droits réservés »

Figure 2 :Figure 3 :Figure 4 :

 234 Figure 2: Anatomical view of the right hemipelvis of a fresh female cadaver. Dissection of the parametrium and paracervix exhibits the duality of the uterine veins and the passage of the ureter. The superficial uterine vein is in front and the deep uterine vein behind.

Figure 5 :

 5 Figure 5: Virtual dissection from the Anatomage® Table of a female left hemipelvis.In green: the urinary system: the bladder, urethra and the ureter. In purple: the uterus and the vagina. In yellow, the inferior hypogastric plexus.

ACKNOWLEDGEMENTS

The authors sincerely thank those who donated their bodies to science so that anatomical research could be performed. Results from such research can potentially increase mankind's overall knowledge that can then improve patient care. Therefore, these donors and their families deserve our highest gratitude (Iwanaga et al., 2021).

The authors would also like to acknowledge Jean-Eric Berton, Julien Landreau and Anatomage® for their technical support.

CA_23850_Figure 1.tiff CA_23850_Figure 2.tiff CA_23850_Figure 3.tiff CA_23850_Figure 4.tiff CA_23850_Figure 5.tiff


Akylas Stratigakos
email: akylas.stratigakos@minesparis.psl.eu

Dennis Van Der Meer
email: dennis.vandermeer@minesparis.psl.eu

Simon Camal
email: simon.camal@minesparis.psl.eu

Georges Kariniotakis
email: georges.kariniotakis@minesparis.psl.eu

End-to-end Learning for Hierarchical Forecasting of Renewable Energy Production with Missing Values

Keywords: data-driven optimization, decision trees, hierarchical forecasting, missing data, renewable energy forecasting

Power systems feature an inherent hierarchical structure. Ensuring the coherency of renewable production forecasts across a hierarchy presents an emerging challenge in energy forecasting. In this context, proposed reconciliation or end-toend learning approaches assume coherent historical observations by construction; this assumption, however, is often violated in practice due to equipment failures. This paper proposes an endto-end learning approach for hierarchical forecasting that directly handles missing values. First, we show that a class of off-the-shelf machine learning models already leads to coherent hierarchical forecasts. Next, we describe a conditional stochastic optimization approach based on prescriptive trees for end-to-end learning with missing values, that fully utilizes the available data. We validate the proposed approach in two case studies comprising 60 wind turbines and 20 photovoltaic parks, respectively. The empirical results show that end-to-end learning outperforms twostep reconciliation approaches and that the proposed solution mitigates the adverse effect of missing data.

I. INTRODUCTION

The increasing penetration of geographically distributed renewable energy sources brings a plethora of challenges in the operational management of power systems. To ensure consistent decision-making, production forecasts at different aggregation levels should adhere to a set of constraints imposed by a hierarchical structure; this presents an increasingly important challenge in energy forecasting.

Traditional approaches for hierarchical forecasting [START_REF] Petropoulos | Forecasting: theory and practice[END_REF], i.e., forecasting a group of time series that satisfy a set of linear aggregation constraints, include the bottom-up and the top-down. The bottom-up approach involves forecasting the bottom-level series and aggregating them; however, it usually performs poorly as the signal-to-noise ratio tends to be lower at the bottom-level series. Further, the top-down approach may introduce forecast bias. Thus, a significant body of research focuses on two-step methods, where each series is modeled independently, with individual (termed base) forecasts being reconciled in a post-processing step. In [START_REF] Wickramasuriya | Optimal forecast reconciliation for hierarchical and grouped time series through trace minimization[END_REF] unbiased base forecasts are reconciled by minimizing the trace of the forecast error covariance matrix. The authors of [START_REF] Van Erven | Game-theoretically optimal reconciliation of contemporaneous hierarchical time series forecasts[END_REF] propose reconciliation by weighted projection of base forecasts. A Parts of this research were carried in the frame of the Smart4RES project (No. 864337), supported by the Horizon 2020 Framework Program. constrained multivariate regression framework is described in [START_REF] Di Modica | Online forecast reconciliation in wind power prediction[END_REF], considering both batch and online learning, applied in wind power forecasting. Moving beyond point forecasts, [START_REF] Taieb | Coherent probabilistic forecasts for hierarchical time series[END_REF] describes a bottom-up approach for coherent probabilistic forecasts, while [START_REF] Yang | Reconciling solar forecasts: Probabilistic forecast reconciliation in a nonparametric framework[END_REF] proposes a block bootstrap method for probabilistic photovoltaic (PV) production forecasts. Recently, end-to-end learning, i.e., training a single model to predict all series in one-shot, has begun to attract attention, as it directly leverages dependencies across series. A deep learning model with an internal projection step is presented in [START_REF] Rangapuram | End-to-end learning of coherent probabilistic forecasts for hierarchical time series[END_REF] for coherent probabilistic forecasts. Lastly, a general framework for end-toend forecasting of predictive quantiles is proposed in [START_REF] Han | Simultaneously reconciled quantile forecasting of hierarchically related time series[END_REF].

A recurring assumption in the literature is that training observations are coherent by construction. In practice, however, missing or erroneous values are commonplace due to communication failures or equipment malfunctions. Two main approaches are identified for dealing with missing values. The first is to ignore observations with missing values (complete case analysis), which is typically applied in the works mentioned above. However, in an end-to-end learning setting, disregarding observations leads to a significant loss of information. Further, if data are not missing at random, bias might be introduced. The second approach involves missing data imputation; in turn, this raises the problem of ensuring that imputed values are coherent. To this end, [START_REF] Liu | Missing value estimation for hierarchical time series: A study of hierarchical web traffic[END_REF] proposes an iterative algorithm to impute missing values while exploiting dependencies across series. In the present work, we examine the case of missing values in the lower levels of the hierarchy, but accurate measurements in the upper levels, which we believe to be of practical interest in power system applications. For example, smart-meters might fail to transmit consumption data at a household level while the respective distribution feeder properly measures aggregated demand.

In short, our contributions are as follows:

• First, we show that a class of non-parametric machine learning models directly derives coherent point and probabilistic hierarchical forecasts, enabling us to use off-theshelf tools in an end-to-end learning setting. • Next, we describe a decision tree algorithm for end-toend forecasts with missing values that does not require imputation and fully utilizes the available training data. • Lastly, we validate the proposed solution in two case studies of day-ahead forecasting considering an aggregation of 60 wind turbines and 20 PV parks, respectively. The rest of this paper is organized as follows. Section II introduces the mathematical background and the proposed methodology. Section III presents the experimental setting and the results. We conclude and provide directions for further research in Section IV.

II. BACKGROUND AND METHODOLOGY

A. Hierarchical forecasting

We examine a group of time series that satisfy a set of linear aggregation constraints, thus forming a hierarchy. Let y it ∈ R denote a realization of the i-th series, associated with features x it ∈ R p , for t = 1, . . . , T . We define y t ∈ R n as the realization of the n series and x t := [x 1t , . . . , x nt ] ∈ R n×p . Following [START_REF] Taieb | Coherent probabilistic forecasts for hierarchical time series[END_REF], let y b t ∈ R n b denote the bottom-level series (i.e., series at the leaf nodes), y a t ∈ R na denote the series formed by aggregation, with n = n b +n a , and S ∈ {0, 1} n×n b denote an aggregation matrix. For a specific hierarchy, observations for all the levels are derived by

y t = Sy b t ⇐⇒ y a t y b t = S a I n b y b t ∀t ∈ [T ], (1)
where S a ∈ {0, 1} na×n b aggregates the bottom-level series y b t , I n b is an n b -size identity matrix, and [T ] := {1, . . . , T }. The notation is illustrated in Fig. 1. Following [START_REF] Rangapuram | End-to-end learning of coherent probabilistic forecasts for hierarchical time series[END_REF], a convenient way to represent (1) is given by

Ay t = 0 ∀t ∈ [T ], (2)
where A = [I na , -S a ] ⊺ and I na is an n a -size identity matrix. Let y t+k denote a set of base point forecasts issued at time t with horizon k, i.e., independent forecasts for each series. The notion of additive coherency (coherency hereafter) is defined as follows.

Definition 1. Forecasts y t+k are said to be coherent if they satisfy A y t+k = 0.

Typically, base forecasts y t+k will not satisfy the coherency constraints, thus a post-processing step is required. Let us further define S := {y | Ay = 0} to be the feasible set that satisfies the linear aggregation constraints. From (1), the following assumption is in place. This assumption is standard in the forecasting literature. In the following, we show that a class of non-parametric machine learning models directly provides coherent point forecasts. First, we state a standard result from convex analysis. Proposition 1. Any convex combination of historical observations y t satisfies the coherency constraints.

Proof. This follows from convexity of S and Assumption 1.

The above holds for additional convex constraints, e.g., nonnegativity of forecasts. A corollary of Proposition 1 is that a class of machine learning models, including, among others, k-nearest neighbors, kernel regression, and decision trees, directly leads to coherent forecasts. These models are based on the idea of local averaging or smoothing of historical observations. For an out-of-sample observation x t+k , we derive a set of non-negative weights ω t (•), with t∈[T ] ω t (x t+k ) = 1, and the respective point forecast y t+k is given by

y a 1 y a 2 y b 1 y b 2 y a 3 y b 3 y b 4 y b 5 S =     1 1 1 1 1 1 1 0 0 0 0 0 1 1 1 I 5    
y t+k = t∈[T ] ω t (x t+k )y t , (3)
i.e., a convex combination of historical observations. From Proposition 1, we see that y t+k are coherent. The practical implication is that off-the-shelf machine learning tools are readily applicable for end-to-end hierarchical forecasting. This result is somewhat trivial; nonetheless, it seems to have escaped the respective forecasting literature.

The above-mentioned models can also be employed for probabilistic hierarchical forecasting. For example, the selected neighbors y t in a k-nearest neighbor model can be treated as (coherent) sample path realizations of the joint predictive density of all the series in the hierarchy. Similarly, one could treat the output of individual trees within an ensemble as realizations of the multivariate predictive density. Therefore, off-the-shelf machine learning tools are also applicable to probabilistic hierarchical forecasting, presenting a computationally cheaper alternative to the internal sampling and projection approach proposed in [START_REF] Rangapuram | End-to-end learning of coherent probabilistic forecasts for hierarchical time series[END_REF] and the bottom-up method described in [START_REF] Taieb | Coherent probabilistic forecasts for hierarchical time series[END_REF].

B. Prescriptive trees for learning with missing values

Next, we examine the case when bottom-level series have missing values due to equipment failures, but aggregated series maintain correct measurements. We assume that missing values are set at 0, therefore Ay t ̸ = 0 and y t ̸ ∈ S. Without loss of generality, the term "missing" refers both to missing and erroneous measurements, as long as these are identified, e.g., by applying an outlier detection filter, and do not propagate through the hierarchy. We examine this problem under a conditional stochastic optimization lens, integrating predictive and prescriptive analytics [START_REF] Bertsimas | From predictive to prescriptive analytics[END_REF], and formulate prescriptive trees for end-to-end hierarchical forecasting. Prescriptive trees [START_REF] Stratigakos | Prescriptive trees for integrated forecasting and optimization applied in trading of renewable energy[END_REF] refer to decision trees that output prescriptions rather than predictions. In our case, and with a slight abuse of terminology, the prescriptions correspond to hierarchical point forecasts, which must satisfy the coherency constraints (and possibly additional ones). Out-of-sample conditional prescriptions are derived via a weighted Sample Average Approximation (SAA) of the original stochastic optimization problem [START_REF] Bertsimas | From predictive to prescriptive analytics[END_REF].

We follow the popular CART method [START_REF] Breiman | Classification and regression trees[END_REF] by recursively partitioning the feature space with locally optimal splits. Mathematically, a node split separates a feature space R ⊆ R n×p at feature j and point s into two disjoint partitions

R = R l ∪ R r , such that R l = {t ∈ [T ] | x tj < s} and R r = {t ∈ [T ] | x tj ≥
s}, with scalar x tj denoting the t-th observation of the j-th feature 1 . Here, we minimize a generic cost function subject to a set of linear aggregation constraints. The problem of finding the locally optimal split is given by

min j,s   min z l ∈S t∈R l (j,s) c(z l ; x t ) + min zr∈S t∈Rr(j,s) c(z r ; x t )   , (4)
with subscripts l, r referring to the left and right child node, z {l,r} ∈ R n being the locally constant decisions (i.e., hierarchical forecasts), which satisfy the coherency constraints S, R {l,r} being index sets, and c(•) being the cost function to be minimized. Thus, the main difference from the CART algorithm is the requirement for predictions to satisfy a set of constraints and the use of a generic, task-based loss function.

In a typical regression setting, c(•) would correspond to the squared ℓ 2 norm. Here, we want to ignore missing values, without disregarding any quality data points. To this end, we employ an indicator matrix Γ ∈ {0, 1} n×T that checks whether historical observations are missing. A single entry γ of Γ is given by

γ it = 1, if y it is missing 0, otherwise ∀i ∈ [n], t ∈ [T ], (5)
and the cost of sample t is given by

c(•) = ∥y t -(1 -γ t ) ⊙ z∥ 2 2 , (6)
with γ t denoting the t-th column vector of Γ and ⊙ denoting the elementwise multiplication. Note that reliable nodes that lead to series without missing values are weighted more heavily in the objective, thus assuming greater importance during training, which we consider to be a desirable property for this application.

As discussed, forecasts are required to satisfy the coherency constraints imposed by S. Note that (4) involves two equality constrained quadratic sub-problems. An analytical solution is derived by solving a system of linear equations obtained from the Karush-Kuhn-Tucker (KKT) optimality conditions (see Appendix A). Other possible constraints, e.g., non-negativity or monotonicity, can be readily included. In this case, a general-purpose convex solver can be called on to evaluate [START_REF] Di Modica | Online forecast reconciliation in wind power prediction[END_REF]. For an out-of-sample observation x t+k point forecasts are derived via a weighted SAA given by

y t+k = arg min z∈S t∈[T ] ω t (x t+k )∥y t -(1 -γ t ) ⊙ z∥ 2 2 . (7)
In general, decision trees are prone to overfitting. Randomization-based ensembles provide a remedy and lead to impressive predictive performance [START_REF] Januschowski | Forecasting with trees[END_REF]; these are readily applicable within the proposed framework, leading to a prescriptive forest. A single tree is fully compiled, with its leaves outputting coherent forecasts. The corresponding weights are:

ω t (x t+k ) = I[R(x t ) = R(x t+k )] |R(x t+k )| , (8)
where R(x t+k ) is the leaf that out-of-sample observation x t+k falls into, | • | the leaf cardinality, and I[•] an indicator function that checks whether training observation x t falls into R(x t+k ). Lastly, for an ensemble of B trees the weights are given by

ω t (x t+k ) = 1 B B b=1 I[R b (x t ) = R b (x t+k )] |R b (x t+k )| . (9)

III. EXPERIMENTAL SETTING AND RESULTS

This section details the experiment and respective results. We examine performance for day-ahead forecasting of wind and PV production, considering only point forecasts of hourly resolution. As this case study involves energy production, we also require forecasts to be non-negative, thus S := {y|Ay = 0, y ≥ 0}; Proposition 1 holds as S remains convex.

1) Data: We use power measurements from 60 wind turbines and 20 PV parks located in mid-west France, with a nominal aggregated capacity of 120 MW and 4 MW, respectively. The available data sets span the period from December 2018 to September 2020 with an hourly resolution. The first 15 months of data are used for training and tuning, with the remaining 5 months used for testing the performance. For the PV data, only measurements recorded at a zenith angle under 85 o are considered, thus the effective sample size is smaller. Table I provides a summary of the two data sets. We examine performance in point forecasting, with a horizon of 12-36 hours ahead, which is typical in electricity market applications. In both cases, a 3-level hierarchy is considered. Wind production data are naturally aggregated at park level (13 wind parks in total); for the PV production data we construct a fictitious hierarchy based on spatial k-means clustering. Fig. 2 provides an overview of the geographical distribution of the power plants.

For feature data we employ a grid of numerical weather predictions (NWPs) obtained from ECMWF, issued daily at 00:00 UTC with a spatial resolution of 0.1 o × 0.1 o . The 5 extracted features are: surface solar radiation downwards, 100m U-and V-wind speeds, 2m temperature, and total cloud cover. Historical production lags are not considered, as these typically do not improve forecasts in the horizon of interest. Further, including historical lags would introduce missing values in x, which is outside of our scope. For the i-th series, respective feature vector x it comprises the NWPs from the closest grid point in terms of Euclidean distance; when forecasting a group of series in one-shot, respective features are concatenated in a single vector.

2) Verification: Individual production data are normalized by the nominal capacity and performance is assessed in terms of normalized root mean squared error (NRSME). Further, following [START_REF] Di Modica | Online forecast reconciliation in wind power prediction[END_REF], the Scaled RMSE (SRMSE) is used by dividing by the total number of child notes. For a set of T forecasts, the SRMSE for the i-th series is defined as:

SRMSE i =   1 T t∈[T ] y it -y it s i 2   1 2 , ( 10
)
where s i the number of child nodes.

3) Benchmark models: We compare the efficacy of a variety of post-processing and end-to-end learning approaches. Note that for post-processing methods, any learning algorithm can be used to generate base forecasts. To allow for a fair comparison, similar base learners are considered in all cases, i.e., randomized tree ensembles based either on the Random Forest [START_REF] Breiman | Random forests[END_REF] or the ExtraTrees [START_REF] Geurts | Extremely randomized trees[END_REF] algorithm. We focus on tree ensembles exclusively, due to their impressive performance [START_REF] Januschowski | Forecasting with trees[END_REF]. The following approaches are examined:

• BASE: Base forecasts with a Random Forest model for each series, without reconciliation. Typically, these will not be coherent. • BASE-BU: Bottom-up reconciliation applied to the base forecasts of the bottom-level series. • BASE-PRJ: Base forecasts post-processed with a Euclidean projection step. The reconciled forecasts are given by arg min y∈S ∥y -

y t+k ∥ 2 , (11)
where y t+k are the base forecasts at time t with horizon k. Alternative reconciliation methods, such as MinT [START_REF] Wickramasuriya | Optimal forecast reconciliation for hierarchical and grouped time series through trace minimization[END_REF] and constrained multivariate least squares [START_REF] Di Modica | Online forecast reconciliation in wind power prediction[END_REF], were also examined. However, as these methods require additional training, results with missing values were not robust and thus are omitted. • EtE: A single Random Forest model predicting the whole hierarchy, to examine the efficacy of end-to-end learning. • EtE-PF: End-to-end learning with prescriptive forests to deal with missing values. In all cases, except for EtE-PF, observations with missing values are disregarded prior to training. For the EtE approach, this means that observation y t is disregarded if at least two of the n series have a missing value, thus this approach deals with the largest loss of information. A grid search is performed to tune the hyperparameters of the Random Forest models. For the EtE-PF, we employ random node splits to speed-up computations, following the ExtraTrees algorithm [START_REF] Geurts | Extremely randomized trees[END_REF], and similarly perform a grid search for tuning.

4) Results: We simulate the effect of missing values due to equipment malfunctions by sampling a subset of bottomlevel nodes and setting a percentage of training observations to zero. We vary both the number of nodes and the percentage of missing values per node and repeat the experiment 5 times to derive aggregate statistics. For simplicity, we assume that missing values occur at the same timestamp for all nodes. Fig. 3 shows the aggregated SRMSE as a function of the number of sampled nodes and the percentage of missing values per node. Overall, the following are observed: i) end-to-end learning (EtE) outperforms post-processing methods but is also more heavily affected by missing values, ii) post-processing methods are robust against missing values, and iii) the proposed EtE-PF combines the best of both worlds. Regarding the wind data set (Fig. 3a), both EtE and EtE-PF show improved accuracy for lower percentage of missing values, with the performance of EtE gradually degrading as the percentage of missing values increases. Conversely, the EtE-PF proves to be robust, consistently outperforming the reconciliation methods. This result persists both for the case of an increased number of malfunctioning nodes and an increased percentage of missing values per node. Further, BASE-PRJ performs, albeit slightly, better than the BASE and BASE-BU, corroborating previous findings on the benefits of post-processing. Similar results are observed for the PV data set (Fig. 3b), which has a smaller sample size. Overall, the relative increase in average SRMSE for EtE from the smallest (5%) to the largest (50%) percentage of missing observations is 0.8% for the wind data set and 0.6% for the PV data set.

By examining the accuracy of end-to-end learning as a function of the number of selected nodes and the percentage of missing observations we observe that the former has a negligible effect in overall performance; this is partly attributed to the design of the the same timestamp across all nodes. On the contrary, the percentage of missing values has a more pronounced effect.

In order to present a comprehensive study, we repeat the above experiment only for EtE-PF with missing values occurring at different timestamps. The results presented in Fig. 4 are similar to the ones achieved before, with the forecast accuracy decreasing only slightly as the number of nodes increases. Thus, we conclude that the proposed EtE-PF successfully mitigates the adverse effects of missing values in the lower levels of the hierarchy. Lastly, we examine performance for each level of the hierarchy. From Table II we observe that in all cases the SRMSE is lower for higher levels of aggregation due to the spatial smoothing effect. The effect is more pronounced for the wind production data, which we partly attribute to the larger number of wind production series examined. For both data sets, the EtE-PF leads to the best performance in the 2nd and 3rd (bottom) level, while BASE-PRJ leads to the best performance for the 1st (top) level, with the results being, generally, statistically significant. Overall, both EtE and EtE-PF consistently improve performance for the bottom-level nodes, highlighting the benefits of exploiting dependencies across time series in an end-to-end learning setting. Note that the results shown in Table II are obtained over all the iterations (for uniform timestamps); the difference between EtE-PF and EtE becomes statistically significant if we only examine adverse scenarios (higher percentage of missing values), as the performance of EtE declines.

IV. CONCLUSIONS

In this work, we examined the problem of forecasting the production of geographically distributed renewable energy sources with missing values under a conditional stochastic optimization lens. First, we showed that a large class of non-parametric machine learning models generates coherent hierarchical forecasts. The main practical implication is that off-the-shelf machine learning tools can be utilized for hierarchical forecasting, presenting an easy way to create benchmarks. Then, we proposed a prescriptive trees algorithm for end-to-end learning with missing values. Performance was evaluated in two case studies of wind and PV production point forecasting on a day-ahead horizon. Overall, end-to-end learning showed improved aggregate performance against twostep reconciliation approaches; for the bottom-level series this improvement was 1.7% and 2.8% for the wind and PV data, respectively. Conversely, reconciliation approaches proved to be more robust against the number of missing values. The proposed solution managed to combine the best of both worlds as it maintained improved performance while also mitigating the adverse effect of missing data for end-to-end learning.

In future work, we plan to study probabilistic hierarchical forecasting with missing values and also examine the effect of network losses on coherency. Another interesting direction is to examine missing or corrupt values in the feature vector.

Assumption 1 .

 1 Historical observations y t are coherent by construction, i.e., y t ∈ S ∀t ∈ [T ].

Fig. 1 :

 1 Fig. 1: Example of a 3-level hierarchy (top) and the respective aggregation matrix (bottom) with n = 8, n a = 3, and n b = 5.

Fig. 2 :

 2 Fig. 2: Overview of the NWP grid points and the wind and PV plants.

Fig. 3 :

 3 Fig. 3: Aggregated SRMSE for the hierarchy as a function of the number of sampled nodes and the percentage of missing values per node over 5 iterations. Bars correspond to one standard deviation.

Fig. 4 :

 4 Fig. 4: Performance of EtE-PF for missing values at different timestamps. The lines indicate the number of malfunctioning nodes, the shaded areas show one standard deviation.

TABLE I :

 I Data set description

		Train / test sample size n b	na	Capacity (MW)
	Wind 13875 / 3671	60	14	120
	PV	6355 / 2322	20	4	4
					39
					:LQG
					1:3
	/DWLWXGH
			/RQJLWXGH

TABLE II :

 II Average SRMSE (± one standard deviation) per hierarchy level. The best-performing model is underlined in bold font. Bold font indicates that a result does not differ from the best-performing model at the 1% level (Welch's t-test).

	Data set Level (# nodes) BASE	BASE-BU	BASE-PRJ	EtE	EtE-PF
		1	0.0969 ± 0.0002	0.0977 ± 0.0001	0.0968 ± 0.0002	0.0972 ± 0.0005	0.0971 ± 0.0003
	Wind	2 (13)	0.1327 ± 0.0132	0.1322 ± 0.0132	0.1326 ± 0.0123	0.1306 ± 0.0128 0.1305 ± 0.0127
		3 (60)	0.1453 ± 0.0144	0.1453 ± 0.0144	0.1452 ± 0.0138	0.1429 ± 0.0141 0.1428 ± 0.0141
		1	0.0761 ± 0.0001	0.0762 ± 0.0001	0.0759 ± 0.0001	0.0765 ± 0.0005	0.0764 ± 0.0002
	PV	2 (3)	0.0812 ± 0.0040 0.0811 ± 0.0041 0.0811 ± 0.0037 0.0808 ± 0.0035 0.0807 ± 0.0035
		3 (20)	0.1008 ± 0.0176	0.1008 ± 0.0176	0.1009 ± 0.0170	0.0980 ± 0.0153	0.0980 ± 0.0153

For brevity of exposition we focus on quantitative features, although it is straightforward to also include categorical features

ACKNOWLEDGMENT

The authors wish to thank HESPUL, SERGIES, and VES-TAS for the provision of renewable production data and ECMWF for Numerical Weather Predictions.

APPENDIX

This section describes analytical solutions for the SAA subproblems in (4). For t ∈ [T ], the respective SAA is given by

For simplicity, the objective is scaled. The KKT optimality conditions for this problem can be written as

where v ∈ R na denotes the dual variables. We write (13) as

where

) is an n-size diagonal matrix whose entries equal the number of non-missing values per series. Hence, we need to solve this set of n + n a linear equations in the n + n a variables. Lastly, note that it is possible for P to become singular; in this case, the least-squares solution of ( 14) can be used.