Laura Rampazzi
email: laura.rampazzi@uninsubria.it

Biagio Rizzo

Chiara Colombo

Claudia Conti

Marco Realini

Ugo Bartolucci

Maria Perla Colombini

Andrea Spiriti

Laura Facchin

The stucco technique of the Magistri Comacini: the case study of Santa Maria dei Ghirli in Campione d'Italia (Como, Italy)

Keywords: stuccoes, Baroque, artistic technique, Magistri Comacini, microscopic analyses, spectrometric analyses, laser ablation

come

Introduction

The study presented in this article investigates the artistic technique used for the Baroque stucco decorations in Santa Maria dei Ghirli, located in Campione d'Italia (Como), a sovereign Italian territory completely surrounded by Swiss territory and facing the lake of Lugano.

The research is part of a wider program of analytical campaigns within the project L'arte dello Stucco nel Parco dei Magistri Comacini -Valorizzazione, conservazione e promozione (The art of stucco in the Magistri Comacini Park -Enhancement, conservation and promotion, referred to as 'Stucco Project' in the text), funded by the European Community Cross-border co-operation III A Program in order to promote the cultural heritage in the borderline area between Italy and Switzerland [START_REF] Authors | Proceedings of the Final Conference on the Interreg IIIA Project 'L'arte dello stucco (The art of stucco)[END_REF]. In particular, the project managed by Insubria University of Como (Italy) and SUPSI University of Lugano (Switzerland) was focused on the stucco works from five churches and some civil buildings spread between the lakes of Como and Lugano, with the aim of comparing the decoration techniques.

The authors of these decorations were the so-called 'Magistri Comacini', known also as 'Artists of the Lakes': painters, sculptors, architects coming from the aforementioned lakes, and working from the Middles Ages to nowadays. The artists used to work all over Italy and Europe for most of the year and return home during wintertime. On these occasions they usually enriched the local heritage, e.g. by making stucco decorations, which are still unknown outside the region and need to be investigated and compared to similar case studies all around Europe.

One of the most outstanding sites decorated by the Magistri Comacini is the church of Santa Maria dei Ghirli, whose plaster decorations date back to the great period of Baroque stucco artworks. In that period, i.e. 17 th and 18 th centuries, stucco had a starring role in architectural decorating, both polychromatic and enriched by gilded finishing, and sometimes made with casts in order to improve the modelling. The mixture of lime or gypsum binder, sometimes combined with sand and rock fragments, gripped the overhanging finishings to iron or wood frames, and the organic additives ensured workability and plasticity [START_REF] Davey | A history of building materials[END_REF]. Many finishing layers could be applied for aesthetic or conservative purposes and a complex stratigraphy is often observed on the external surface. The Scagliola technique, very similar to stucco technology, was perfected in Emilia Romagna, Tuscany and Lombardy (Intelvi Valley) in the same period. Altarpieces, sculptures, and panels were produced mixing gypsum, organic additives and pigments in order to imitate precious natural stones such as marbles.

New light on the stucco technique used in Santa Maria dei Ghirli was shed by the exploitation of an analytical protocol set up in a previous work [START_REF] Rampazzi | The stuccoes of St. Lorenzo in Laino (Como, Italy): the materials and the techniques employed by the "Magistri Comacini[END_REF] and by the discussion of the results with art historians, in order to be sure that they make sense if compared to the cultural and artistic scenario of that time and to the artist's fingerprint. This paper describes an interdisciplinary and complementary team approach, which resulted in enhancing our knowledge of the decoration from the point of view of its conservation state, composition, technique of preparation, history, iconography, philology, chronology, and relationships to contemporary artists.

In particular, the analytical protocol [START_REF] Rampazzi | The stuccoes of St. Lorenzo in Laino (Como, Italy): the materials and the techniques employed by the "Magistri Comacini[END_REF]) made use of various analytical techniques, complementary in decrypting the layer arrangement and the composition of the stucco decorations: Fourier transform infrared spectroscopy (FTIR), X-ray diffraction (XRD), scanning electron microscope equipped with energy dispersive spectroscopy (SEM-EDS), gas chromatography-mass spectrometry (GC-MS) and laser ablation mass spectrometry equipped with plasma source (LA-ICP-MS). The aforementioned techniques are some of the most popularly used methods. FTIR identifies the inorganic and organic compounds and XRD determines the crystalline phases of the minerals [START_REF] Montana | The 'recipe' of the stucco sculptures of Giacomo Serpotta[END_REF][START_REF] Cavallo | The emigration of prealpine lake masters: a study of a 17th-century mantelpiece attributed to the Carlone family of Rovio[END_REF][START_REF] Tulliani | Study of the degradtion causes affecting stucco sculptures from the Valentino Castle in Turin[END_REF]Cardell-Fernandez 2006;[START_REF] Galván-Ruiz | Chemical and physical characterizationof stuccos from a Mexican colonial building: el Museo del Calendario of Queretaro[END_REF][START_REF] Mazzocchin | Analysis of roman wall paintings from the thermae o f 'Iulia Concordia[END_REF][START_REF] Sansonetti | Colored grounds of gilt stucco surfaces as analyzed by a combined microscopic, spectroscopic and elemental analytical approach[END_REF]); SEM-EDS investigates the morphology and maps the element along the stratigraphy [START_REF] Montana | The 'recipe' of the stucco sculptures of Giacomo Serpotta[END_REF][START_REF] Galván-Ruiz | Chemical and physical characterizationof stuccos from a Mexican colonial building: el Museo del Calendario of Queretaro[END_REF][START_REF] Sansonetti | Colored grounds of gilt stucco surfaces as analyzed by a combined microscopic, spectroscopic and elemental analytical approach[END_REF]; GC-MS allows the determination of the organic content [START_REF] Wouters | Baroque stucco marble decorations. II. Composition and degradation of the organic materials in historical samples and implications for their conservation[END_REF] and LA-ICP-MS traces the element content on the surface [START_REF] Giussani | Role of Laser Ablation-Inductively Coupled Plasma-Mass Spectrometry in Cultural Heritage research: a review[END_REF]. In particular, the coupling of LA-ICP-MS and SEM proved effective in distinguishing layers very thin and similar from the compositional point of view, thanks to the high sensitivity for trace elements analysis of the former and to the high resolution of the latter [START_REF] Rampazzi | The stuccoes of St. Lorenzo in Laino (Como, Italy): the materials and the techniques employed by the "Magistri Comacini[END_REF].

Historical and artistic background

The cycle of the Santa Maria dei Ghirli church is an interesting case of stucco decoration in early Baroque style, still influenced by mannerist accents. The building is placed in Campione d'Italia, a small Italian enclave in Switzerland that actively contributed toward the artistic development of many European cities since medieval times through the work of the Campione Masters, part of the wider group of the 'Artists of the Lakes' (A. Spiriti in Various Authors, 2006, pages 27-32). The main nucleus of Santa Maria dei Ghirli church was probably built in the 7-8 th century, and between 777 and 874 it became the property of the monastery of Sant' Ambrogio in Milan and under its jurisdiction. The growing devotion to the Virgin Mary (the church is dedicated to the Annunciation) favoured enlargements and remakings of the building over the centuries, particularly in the second half of the 14 th and early 15 th century [START_REF] Colombo | Storia del monumento, in 'Il Santuario di Santa Maria dei Ghirli in Campione d'Italia[END_REF]).

In the first decades of the 17 th century important architectural and decorative modifications gave to the church the present aspect. The project was planned by Isidoro Bianchi (Campione d'Italia, 1581-1662), an artist who used to work as painter, sculptor and engineer in Prague, Lombardy and especially as court painter for the Savoy dynasty in Turin where he was the leader of a team which involved firstly Antonio and Pietro Castelli and later his sons, Pompeo and Francesco [START_REF] Grabach | Allgemeines Künstler-Lexicon[END_REF][START_REF] Dell'omo | Da Praga a Torino: gli inizi di Isidoro Bianchi. Nuove ipotesi e qualche novità dai documenti[END_REF]. According to the dates painted and carved on the walls, the work started at least in 1610-12 and lasted until 1634 [START_REF] Angelis | Affreschi e stucchi del XVII secolo, in GIl Santuario di Santa Maria dei Ghirli in Campione d'Italia[END_REF]. Isidoro Bianchi planned the new lantern and dome connected with the third span of the single nave and covered by a barrel vault, partially remodelled the presbytery area, and coordinated the artist's team that dealt with the interior decoration, i.e. frescoes and stuccoes.

The Ghirli cycle had been only considered in a wider range of studies on Lombardy religious building and on Isidoro Bianchi himself [START_REF] Colombo | Storia del monumento, in 'Il Santuario di Santa Maria dei Ghirli in Campione d'Italia[END_REF][START_REF] Angelis | Affreschi e stucchi del XVII secolo, in GIl Santuario di Santa Maria dei Ghirli in Campione d'Italia[END_REF], but it had never been the object of a specific essay until the 2 Experimental Some meaningful decorations (Figure 1) were sampled according to art historians' suggestions (Table 1). In particular, the historical and artistic background influenced the samples' selection on two principles. The former suggested the localization in those areas where Bianchi's work was documented, i.e. the presbytery, the vault, the high altar and at the highest part of the nave vault. The second principle suggested sampling different typologies of stucco patterns, i.e. the decorative parts (e.g. the scrolls, normally executed by the master's team with use of moulds) and the altar angels'

statues with an in the round decoration, certainly made by the lead artist Isidoro Bianchi.

Results and discussion

Information about the stratigraphy was obtained from observations carried out on thin and polished cross sections by optical and scanning electron microscopy, while EDS, XRD, FTIR, GC-MS and LA-ICP-MS techniques clarified the composition, in terms of major, minor and trace elements.

The stucco was generally composed of at least three layers (UNI10924, 2001):

-the inner one, i.e. the bulk, composed by a whitish mortar;

-the middle one, made by a more plastic white mortar showing a thickness around 8 mm;

-the outer part, made mainly by a binder fraction, which was characterised by a white-greyish colour with a thickness around 20-40 microns.

The main components of the bulk was lime binder. Clots were observed within the binder fraction, due to the coagulation of particles during the lime slaking process. The phenomenon was probably caused by the presence inside the slaked lime of a fraction of lime already carbonated before use [START_REF] Bugini | The presence of white lumps in ancient mortars[END_REF]. The aggregate showed a low sorting (with size ranging from few microns to 1-2 mm) and a high packing, and was mainly constituted by fragments of silicate metamorphic rocks, many crystals of silicate minerals and calcite, and by iron oxides. From the mineralogical point of view no significant differences were detected in the samples, as probably coming from the same raw materials sources.

As far as the middle portion of the stuccoes is concerned, the open total porosity values ranged from 25% to 28%, with a prevalent porosimetric distribution in the 0.25 and 0.05 m pores class. The average thickness was about 8 mm. The binder was composed by magnesic lime, as shown by microscopic observations. Calcite and magnesite (MgCO 3 ) signals were determined in FTIR and XRD spectra as well. Samples 1, 3 and 4 also contained hydromagnesite

[Mg 5 (CO 3 ) 4 (OH) 2 4(H 2 O)], as suggested by the peculiar FTIR absorbance peaks at 3649, 3512 and 3447 cm -1 , and 1428, 1423 cm -1 , due to the presence of two kinds of carbonate ions in the molecule [START_REF] Bruni | Characterization of ancient magnesian mortars coming from northern Italy[END_REF]) (Figure 2a). XRD analyses confirmed the presence of hydromagnesite in sample 3 and sample 4. The presence of hydromagnesite is due to the use of dolomite rocks as binder and/or as aggregate (Atzeni 1996). Moreover, according to the literature its crystallization among sand grains usually induces better mechanical strengths to mortars [START_REF] Caceres | Thermal decomposition of dolomite and the extraction of its constituents[END_REF][START_REF] Dheilly | Hydromagnesite development in magnesian lime mortars[END_REF]. Gypsum was determined by thin section observations, XRD and FTIR analyses only in sample 2 and, in traces, in sample 1. Clots of lime and gypsum were observed in samples 2 and 3, as clearly shown by Figure 3a. As regards the aggregate of the middle portion, it showed a high sorting (with average size of 500m) and a low packing, and was mainly constituted of fragments of carbonaceous rocks (both sedimentary and metamorphic) and of low amounts of metamorphic rocks, except for sample 4, that contained mainly carbonatic fragments. XRD investigation ascertained the presence in sample 1 and 2 of phyllosilicates, that may be ascribed to the metamorphic rocks fraction. EDX maps of element distribution of calcium, magnesium, silicon and sulphur showed how the presence of Si in the aggregate was very low, while the gypsum was present both in the binder and in the clots (Figure 3b).

Of significant interest from the investigation of the finishing layers is that the outer part is very simple and only in two samples composed by one or two layers (Figure 4). It is surprising, since according to the technology commonly adopted by Magistri Comacini the stucco decorations were completed with a number of finishing layers [START_REF] Cavallo | The emigration of prealpine lake masters: a study of a 17th-century mantelpiece attributed to the Carlone family of Rovio[END_REF][START_REF] Rampazzi | The stuccoes of St. Lorenzo in Laino (Como, Italy): the materials and the techniques employed by the "Magistri Comacini[END_REF]). The external surfaces of stucco portions of samples 3 and 4 appeared covered only by deposit particles, composed by silicates, iron and calcium oxalate in the former and phosphates in the latter. LA-ICP-MS profile confirmed the enrichment of iron on the external surface of sample 3. Sample 1 presents only one finishing layer, made of magnesic lime, gypsum and calcium oxalate. The composition of sample 2 is peculiar and differs from the other ones. It presents two finishing layers (Figure 5). The external one is composed only by magnesic lime, while the internal one is constituted by lime and gypsum clots. SEM-EDX analyses pointed out the presence of barium, sulphur and zinc, suggesting the presence of barium white (BaSO 4 )

and zinc white (ZnO). Both the pigments date to the 19 th century [START_REF] Feller | Barium Sulphate -Natural and Synthetic in Artists' pigments[END_REF][START_REF] Kuhn | Zinc White in Artists' pigments[END_REF]), thus indicating an ante quem period for the application of a finishing layer, probably on occasion of a conservation work. A restoration of the whole building, including stucco decorations, was actually documented between 1910 and 1911 by the parish priest's report during the periodical pastoral visitation [START_REF] Storico | Magnesian Limes. Experimental contribution to interpreting historical data, Science and Technology for Cultural Heritage[END_REF]. Past surface treatments usually consisted of maintenance application of opaque or glossy coatings to make uniform and protect the decayed surface, or of making the appearance of the ancient decorations and brand new ones uniform [START_REF] Zanardi | La superbia della bianchezza del marmo in Proceedings of the 1st International Symposium The oxalate films. Origin and significance in the conservation of works of art[END_REF]). As far as the characterisation of the organic components of the finishing and stucco layers is concerned, Table 2 reports the amino acid percentage contents (w/w%) obtained by quantitative GC-MS analyses. The results showed that:

-only samples 1 (0.4 ppm both for stucco layer and finishing layer) and 4 (0.4 ppm for stucco layer, 1.6 ppm for external surface) contained amounts of a proteinaceous material above the quantification limit (LOQ =0.05 ppm);

-samples 2 and 3 had amino acid values close to the detection limits (LOD= 0.01 ppm);

-the absence of hydroxyproline, which rules out the presence of animal glue;

-a quite high amount of glycine (15-20%), which can be attributed to the effect of biological agents such as fungi or bacteria: the metabolism product of many bacteria is indeed glycine [START_REF] Colombini | Analytical strategies for characterizing organic paint media using gas chromatography/mass spectrometry[END_REF]. Since the principal component analysis (PCA) of amino acid contents is a very powerful technique for pattern recognition (Colombini 2004), it was applied to the data of Table 2 for samples 1 and 4 which presented a proteinaceous content higher than the LOQ. Figure 6 reports the score plot and the loading plot of PCA applied to the amino acid percentage contents: the first two principal components account for 74.1% of the total variance. The samples form a cluster close to the composition of egg. Basically, the difference in the PC1 values are mainly due to the high percentage of glycine which causes a shift towards the animal glue cluster. Moreover, the effect of the ageing acts on the values of PC2 which are slightly higher in respect to the reference values. Notwithstanding these, the results suggest the presence of egg both in stucco and finishing layers. Otherwise, organic compounds could have been used as polishing surface treatments. The hypothesis seemed to be confirmed by the presence of egg in samples 1 and 4 as well as by the occurrence of apatite in sample 4, as suggested by the FTIR peaks at 1092 and 1034 cm -1 and by SEM-EDX analyses (Figure 7). Apatite may in fact come by the mineralization of organic compounds treatments containing phospho-compounds [START_REF] Rampazzi | Analytical investigation of calcium oxalate films on marble monuments[END_REF]) and as a consequence could have been originated in the same way as calcium oxalate. In samples 1 and 4 a lead enrichment in the external portion of the finishing layers was identified by LA-ICP-MS analyses (Figure 8). The high sensitivity of the technique is effective for the analysis of trace element profiles not detectable by SEM-EDX investigations, thus completing the element maps along the stratigraphy. The polished cross sections of the samples were analyzed by SEM-EDS, in order to determine the major and minor elements and to define into details the stratigraphy, in particular the number and thickness of the layers. In the next step, this map guided the laser along the stratigraphy, linking the distribution of elements to the single layers. The presence of lead suggested the use of lead moulds which "contaminated" the surfaces during the set, leaving trace concentrations of lead. The use of moulds for repeated decorative elements is described in art treatises.

Sample

Figure 8 LA-ICP-MS analysis of sample 1: lead content variation along the stratigraphy.

Conclusions

The analytical campaign on the stucco artworks of S. Maria dei Ghirli church allowed the construction of an hypothesis on Isidoro Bianchi's decoration technique. The comparison between the compositional data of artworks and their semantic and historical context proved effective for the clarification of chronological, and stylistic, problems and the knowledge of materials. In particular, the stratigraphy was able to shed a light on some chronological problems.

According to literary sources, Isidoro Bianchi's team decoration campaign took place in a single phase from 1620 to 1623/25 on the occasion of the 1625 Jubilee year. This is confirmed by the similar composition and stratigraphy of the stucco decorations released in the altar and in the nave (samples 1, 3 and 4). The sample coming from the presbytery vault (sample 2) gave different results for the presence of gypsum in the binder fraction, probably in order to emphasize an illusionistic use of stucco, typical of Baroque style, as imitation of more precious marbles, thanks to the brightness effects of the surface.

The presence of various collaborators in Santa Maria dei Ghirli site, according to 'Artists of the Lakes' tradition, was probably due to a specific request for a quick execution of the interior stucco decoration and to the hurry in accomplishing the campaign before the Jubilee event. This hypothesis is corroborated by the analytical results. The presence of organic compounds enhanced in fact the workability. Moreover, the stratigraphy investigation showed the rather coarse quality of stucco elements and way of shaping, since the finishing layers were surprisingly absent or plain.

The effort to give a bright appearance to the decoration surface suggests that Isidoro Bianchi and his team probably intended to conclude the campaign in a more suitable moment, with much time available, or under more propitious economic conditions. The lead enrichment in the finishing layers of one sample suggested the use of moulds, once again to make the decorating quicker, especially for decorative parts such as cartouches, vegetable elements (i.g. leaves), scrolls.

Notwithstanding the uniformity of chemical and mineralogical composition in 3 of the 4 samples analyzed, some differences in patterns and shaping, together with the documentary evidences, allowed us to distinguish four different artists, probably including Isidoro Bianchi. The identity of Isidoro Bianchi's collaborators is still unknown, except for his sons. It is possible to suggest the name of Antonio and Pietro Castelli from Melide, Isidoro's partners when he worked in Turin. The determination of peculiar pigments dating back to the 19 th century is probably connected to the restoration works between 1910 and 1911. Although not described in the archive documents, the stucco restoration could have been carried out in order to refresh the aged surfaces only where necessary, possibly with a thin layer of plaster as a few traces were found. The most seriously damaged area of the stucco decorations was certainly that of the presbytery vault, from which sample 2 was taken, because of a documented deep crack.

  research carried out within the 'Stucco Project' (various authors 2006). The archival investigation was focused especially on notarial deeds (such as contracts, wills and proxies) kept in the State Archive of Milan, in order to better reconstruct the decoration cycle's chronology and the artists that modelled it. The research moved to the text written by Roberto Rusca (1566-1629/after), Ambrosian monk and Vicar of Campione (1620-1623), most probably the patron of the whole enterprise. His chronicle, published in two different editions (Bergamo 1625 and Piacenza 1629), recorded the execution of most of the stucco decoration during the years 1620-21-22 (Rusca 1625).The other two stucco campaigns were probably carried out by still unknown artists. The former took place in 1686, as suggested by the date found during the last conservation work in winter-spring 2010, involving the frames in the presbytery side walls. The latter dated back to the 3 rd -4 th decade of the 18 th century, as confirmed by archive documents, contemporary to the completion of the façade on Isidoro Bianchi's project. The campaign concerned the presbytery pilaster strips which were decorated with the brand new Rococo sensibility spread throughout Europe by the 'Artists of the Lakes' during the first decades of the 18 th century[START_REF] Spiriti | La nascita del Rococò problemi e certezze, 'Il Rococò: nascita di un linguaggio artistico[END_REF].

Figure 1

 1 Figure 1 Clockwise from upper left: Picture of the church of Santa Maria dei Ghirli, church plan with sample location, photographs showing the stucco decorations from which sample 3 and sample 2 were taken.

Figure 2 FTIR

 2 Figure 2 FTIR transmittance spectra of the stucco layer of sample 4 (a), showing the co-presence of calcite (874, 713 cm -1 ), magnesite (744 cm -1 ) and hydromagnesite (3649, 3512, 3447, 1428 and 1423 cm -1 ), and of the finishing layer of sample 1 (b), showing the presence of weddellite (1647 and 1324 cm -1 ).

Figure 3 .

 3 Figure 3. Thin cross section (a) of sample 2 (bar 200 m), showing the presence of lime and gypsum clots (Nicol 45°), and polished cross section (b) of sample 2 analysed by SEM-EDS, backscattered image (bar 100 m) and EDX distribution maps of Ca, Mg, Si and S of the same area (clockwise from upper left).

Figure 4 .

 4 Figure 4. Polished cross section and stratigraphic scheme of sample 1 (bar 100 m) showing the presence of just one finishing layer.

Figure 5 .

 5 Figure 5. Polished cross section of sample 2 analysed by SEM-EDS. Backscattered image of polished cross section (bar 100 m) and EDX distribution maps of Ca, Mg, Zn, Ba and S of the same area (clockwise from upper left).

Figure 6

 6 Figure 6 PCA score plot (left) and loading plot (right) of amino acid percentage contents of samples from St. Maria dei Ghirli.

Figure 7

 7 Figure 7 Polished cross section of sample 4 analysed by SEM-EDS. Backscattered image of polished cross section (left, bar 50m), EDX distribution maps of P and Ca of the same area (right).

Table 1

 1 Sampling and description of the microfragments of stucco.

Table 2 :

 2 Amino acid percentage contents of the samples and the content (ppm) of proteinaceous material. [ala (alanine), gly (glycine), val (valine), leu (leucine), ile (isoleucine), ser (serine), pro (prolinae), phe (phenilalanine), asp (aspartic acid), glu (glutammic acid), hyp (hydroxyproline), LOD=limit of detection].

		ala	gly	val	leu	ile	ser	pro	phe	asp	glu	hyp	ppm
	1 finishing layer	16.5	19.8	10.6	14.1	7.7	4.9	5.3	5.5	6.7	8.9	0.0	0.4
	1 stucco layer	12.8	16.4	13.6	16.6	5.7	3.3	4.4	2.2	12.9	12.1	0.0	0.4
	2 finishing layer	25.2	24.3	10.1	13.3	16.5	4.3	2.2	2.7	0.8	0.6	0.0	<LOD
	2 stucco layer	23.6	24.3	11.7	13.0	15.8	5.3	0.8	2.3	2.6	0.6	0.0	<LOD
	3 finishing layer	22.0	26.0	12.0	15.2	15.3	3.6	1.2	2.0	2.2	0.5	0.0	<LOD
	3 stucco layer	21.7	26.3	14.2	13.0	18.3	1.6	1.2	1.0	1.2	1.0	0.0	<LOD
	4 external surface	10.2	17.4	12.6	17.3	7.7	4.8	4.8	3.5	12.3	9.3	0.0	1.6
	4 stucco layer	14.5	14.6	14.5	16.8	5.2	3.2	6.8	1.6	15.3	7.6	0.0	0.4

Acknowledgments

The authors gratefully acknowledge the financial support provided by the EU Interreg IIIA Program and D. Di Tommaso and Don E. Mosca for the church plan.


Gratien Bonvin

Sophie Demassey

Claude Le Pape

Nadia Maïzi

Vincent Mazauric

Alfredo Samperio

A Convex Mathematical Program for Pump Scheduling in a Class of Branched Water Networks

Keywords: pump scheduling, branched water networks, demand-response, Mixed-Integer Non-Linear Programming, convex relaxation

We address the day-ahead pump scheduling problem for a class of branched water networks with one pumping station raising water to tanks at different places and levels. This common class is representative of rural drinking water distribution networks, though not exclusively so. Many sophisticated heuristic algorithms have been designed to tackle the challenging general problem. By focusing on a class of networks, we show that a pure model-based approach relying on a tractable mathematical program is pertinent for real-size applications. The practical advantages of this approach are that it produces optimal or near-optimal solutions with performance guarantees in near real-time, and that it is replicable without algorithmic development. We apply the approach to a real drinking water supply system and compare it to the current operational strategy based on historical data. An extensive empirical analysis assesses the financial and practical benefits: (1) it achieves significant savings in terms of operation costs and energy consumption, (2) its robustness to dynamic pricing means that demand-response can be efficiently implemented in this type of energy-intensive utility.

Introduction

In typical drinking water distribution networks, water is collected, treated, and pressurized at pumping stations, then raised to elevated tanks and water towers, and finally supplied by gravity to local customers [START_REF] Swamee | Design of water supply pipe networks[END_REF].

In developed countries, water utilities are among the largest consumers of energy in municipalities, and nearly all of the electricity purchased is used for pumping [START_REF] Denig-Chakroff | Reducing electricity used for water production: Questions state commissions should ask regulated utilities[END_REF]. The motives for installing computer control systems in such water utilities have been documented for decades [START_REF] Shamir | Computer applications for real-time operation of water distribution systems[END_REF]: high operating costs, system complexity, rising demand, need to capitalize on staff experience, etc. The conjugated growth of electricity prices and environmental awareness has spurred new interest in the optimal control problem, also known as the Pump Scheduling Problem, because of its ability to address energy efficiency mechanisms such as demand-response [START_REF] Ashok | Load-management applications for the industrial sector[END_REF][START_REF] Ghaddar | A lagrangian decomposition approach for the pump scheduling problem in water networks[END_REF]. Because they are energy intensive and thanks to their buffering abilities, drinking water distribution networks are indeed remarkable instruments for implementing demand-response programs via load-shifting: by storing water and energy in elevated tanks, $ This paper was presented at the 7th International Conference on Applied Energy (ICAE2015), March 28-31, 2015, Abu Dhabi, UAE (Original paper title: "A heuristic approach to the water networks pumping scheduling issue" and Paper No.: 597).

* Corresponding author Email address: gratien.bonvin@mines-paristech.fr (Gratien Bonvin) demand can be anticipated by pumping the adequate amount of water when electricity is cheapest. Hence, given reliable water consumption forecasts and dynamic electricity prices, significant energy and financial savings can be 25 achieved by shifting the pump operations to off-peak [START_REF] Ghaddar | A lagrangian decomposition approach for the pump scheduling problem in water networks[END_REF][START_REF] Van Staden | A model predictive control strategy for load shifting in a water pumping scheme with maximum demand charges[END_REF].

Due to the problem's complexity, which is even greater when prices are highly dynamic, only automatized control programs have the computational capacity to determine a least-cost pump schedule in a general real-size utility. This 30 issue has given rise to abundant literature on sophisticated heuristic algorithms (see e.g. the surveys in [START_REF] Shamir | Computer applications for real-time operation of water distribution systems[END_REF][START_REF] Ormsbee | Optimal control of water supply pumping systems[END_REF][START_REF] Nicklow | The ASCE task committee on evolutionary computation in environmental and water resources engineering. State of the art of genetic algorithms and beyond in water resources planning and management[END_REF][START_REF] Price | Iterative linearization scheme for convex nonlinear equations: Application to optimal operation of water distribution systems[END_REF][START_REF] Ghaddar | A lagrangian decomposition approach for the pump scheduling problem in water networks[END_REF]). However, many small to medium-sized water supply systems, such as the French rural distribution network we examine in this paper, still rely on human operators to 35 manage pumps in near-real time and expect them to do so in a cost-effective way [START_REF] Feldman | Aspects of energy efficiency in water supply systems[END_REF]. One probable reason is that these practical solutions, including commercial products [START_REF] Bunn | Application of an expert system to control treated water distribution[END_REF][START_REF] Engineering | Optimising water resource and supply systems[END_REF]13], require a minimum of algorithmic development and expertise to be replicated. Furthermore, these 40 heuristics do not guarantee optimality and, thus, cannot be used to assess the maximum potential financial benefit.

By focusing on a specific, still widespread, class of networks, we show in this paper that a pure model-based approach is possible to tackle real applications, with several practical advantages. First, this approach is easy to replicate on any network in the class, since only the mathematical model's input data need to be changed; the resulting model can then be handled by any certified mathematical programming solver. This feature is convenient for net-work operators who wish to quickly prototype the optimal control of their utilities to estimate the benefit. Second, the model is tractable and can be quickly solved. Our control approach can then even be considered to operate in near real-time. Third, the approach is exact, i.e. it can guarantee the maximum possible financial benefit. Finally, since it is robust to highly-dynamic electricity prices, it is suitable for implementing demand-response.

The networks considered in this paper have the following characteristics: (1) a branched layout (i.e. no loop, no bi-directional pipes), [START_REF] Denig-Chakroff | Reducing electricity used for water production: Questions state commissions should ask regulated utilities[END_REF] several pumps collocated in one pumping station at the source node, and (3) elevated tanks equipped with flow-control valves at the destination nodes. While such a topology can be encountered in industrial usages or in irrigation systems [START_REF] Pulido-Calvo | Selection and operation of pumping stations of water distribution systems[END_REF], it is particularly common in rural drinking water distribution networks [START_REF] Swamee | Design of water supply pipe networks[END_REF]: in the upstream pumping network, water is treated and pumped near the source to be raised and dispatched to water towers where valves compensate for the elevation differences. The downstream gravity-fed network -from the towers to the consumers -is independent since the pressures in both networks are not correlated. To our knowledge, this class of networks has not been specifically investigated in prior studies. Notably, benchmark problems and real applications in literature tend to concern urban networks with looped configurations.

Our approach is based on a traditional formulation of the pump scheduling problem as a Mixed-Integer Non-Linear program (MINLP). This non-convex model has not been solved in reasonable time in general cases even involving small networks [START_REF] D'ambrosio | Mathematical programming techniques in water network optimization[END_REF]. We thus consider a tractable relaxation by dropping non-convexities, precisely, the strict equality in the potential-flow coupling constraints. For the considered class of networks, we show how to convert a solution of this relaxation into a feasible solution with a limited cost increase. The increase is even null (i.e. the solution remains optimal) when the characteristic power functions of the pumps share the same slope.

The study was conducted by an industrial application on the FRD network, a medium-sized drinking water distribution system (6 pumps, 16 tanks, 49 nodes, 53 arcs) located in a rural zone of France and supplying about 3 million m 3 of water per year. Experiments on one year of historical data show the practical benefits of the approach compared to the in-house strategy currently employed in 95 the utility: pumping plans for the day-ahead are generated in less than one minute and their feasibility is assessed on the EPANET simulator. On average, over one year, their cost is guaranteed to be less than 3% of the optimum, and they are 17% cheaper than the plans generated heuristically by the current strategy. In addition, the daily plans involve lower discharge pressures (22% lower on average) and thus fewer leakages. Incidentally, they operate the pumps with a 10% lower mean efficiency, which indicates that the current design of the station is sub-optimal. Finally, the current strategy is driven by night/day electricity tariffs and is not compatible with an intra-day price variability, unlike our approach.

The paper is organized as follows: Section 2 provides a selected review of literature on the Pump Scheduling 110 Problem. Section 3 presents the traditional MINLP formulation and Section 4 the convex relaxation for the considered class of networks. Section 5 provides an extended empirical analysis on the FRD network. Section 6 gives conclusions and perspectives.

Literature review

Pump Scheduling involves planning the pumping operations over a future period, typically the day ahead, in order to minimize energy costs, considering the static physical characteristics of the network elements (pipes, 120 pumps, valves, tanks) and forecasts of the electricity price and water demand. An accurate mathematical model for this problem is hard to optimize as it involves both nonlinearities (the potential-flow relations) and combinatorial aspects (i.e. when active elements, pumps and valves can 125 be switched on/off) [START_REF] D'ambrosio | Mathematical programming techniques in water network optimization[END_REF][START_REF] Zhang | Minimizing pump energy in a wastewater processing plant[END_REF].

This problem has been the focus of significant attention over the last 40 years. Several survey papers [START_REF] Shamir | Computer applications for real-time operation of water distribution systems[END_REF][START_REF] Ormsbee | Optimal control of water supply pumping systems[END_REF][START_REF] Price | Iterative linearization scheme for convex nonlinear equations: Application to optimal operation of water distribution systems[END_REF] have reviewed the state-of-the-art. Altogether, they cite more than 80 references that we summarize here: earlier 130 works mostly applied dynamic programming or continuous linear/non-linear programming techniques, often combined with hydraulic simulators, to handle the complex physics out of the optimization models. These time-consuming methods either consider basic networks (one pump, 135 one tank) or oversimplify the problem.

Later on, many metaheuristic procedures were developed, including a number of evolutionary algorithms reviewed in the comprehensive survey [START_REF] Nicklow | The ASCE task committee on evolutionary computation in environmental and water resources engineering. State of the art of genetic algorithms and beyond in water resources planning and management[END_REF] and, recently, less traditional procedures such as particle swarms [START_REF] Zhang | Minimizing pump energy in a wastewater processing plant[END_REF] or mat-140 heuristics, i.e. mathematical programming-based heuristics [START_REF] Puleo | Multi-stage linear programming optimization for pump scheduling[END_REF]. Several works report good results on medium-size benchmark instances or industrial applications: they are computationally efficient and can handle complex facets of the problem (e.g. multiple objectives [START_REF] Savić | Multiobjective genetic algorithms for pump scheduling in water supply[END_REF] or side con-145 straints). On balance, they do not guarantee global optimality and require sophisticated algorithm developments.

To tackle these drawbacks, model-based approaches rely on mathematical formulations to be solved using any suitable general-purpose solvers. Recent advances in Mixed-150 Integer Non-Linear Program (MINLP) solvers have renewed the interest in these approaches, as reported in a very recent survey [START_REF] D'ambrosio | Mathematical programming techniques in water network optimization[END_REF]. However, the global (even local) solvers that are currently available do not yet scale up on the hard MINLP formulation of the pump scheduling 155 problem. Hence, the approaches based on mathematical programming in current literature rely on relaxed models. Three kinds of relaxations are considered: the discrete decision components are ignored [START_REF] Burgschweiger | Optimization models for operative planning in drinking water networks[END_REF][START_REF] Skworcow | Pump schedules optimisation with pressure aspects in complex large-scale water distribution systems[END_REF]; the nonconvexities are approximated using piecewise linear func-dropped [START_REF] Gleixner | Towards glob-860 ally optimal operation of water supply networks, Numerical Algebra[END_REF]. The resulting programs (a continuous nonconvex program, a mixed integer linear program, or several small independent MINLPs) are more tractable, although they remain NP-hard. A trade-off is thus required between the strength of the approximation (e.g. the number of segments in the piecewise linear fit) and the expected computation time. Furthermore, these approaches rely on a heuristic algorithm to recover the feasibility of the relaxed solutions in a second phase. As a direct evolution of relaxation techniques, two recent works apply decomposition techniques (lagrangian relaxation [START_REF] Ghaddar | A lagrangian decomposition approach for the pump scheduling problem in water networks[END_REF] and Benders decomposition [START_REF] Naoum-Sawaya | Simulationoptimization approaches for water pump scheduling and pipe replacement problems[END_REF]): to get an optimal solution, the relaxation is solved iteratively driven by a master program. Again, as the convergence is too slow, the method is truncated and coupled with Local Search [START_REF] Ghaddar | A lagrangian decomposition approach for the pump scheduling problem in water networks[END_REF]. Hence, the mathematical programming-based approaches that were published for the Pump Scheduling Problem are not pure modelbased approaches: they are generic and some are applicable to large networks, but they generally do not provide optimality proofs, many are time-consuming, and all require more or less experimented skills to be implemented properly. In comparison, we focus in this article on a less general (yet widespread) class of networks but, precisely, we show how to exploit these specificities to derive a fast 185 exact replicable pure model-based approach.

MINLP for Pump Scheduling

The Mixed-Integer NonLinear Programs designed for the general Pump Scheduling Problem [START_REF] Burgschweiger | Optimization models for operative planning in drinking water networks[END_REF][START_REF] Gleixner | Towards glob-860 ally optimal operation of water supply networks, Numerical Algebra[END_REF][START_REF] D'ambrosio | Mathematical programming techniques in water network optimization[END_REF][START_REF] Ghaddar | A lagrangian decomposition approach for the pump scheduling problem in water networks[END_REF] share the same structure. For the class of networks considered 190 in this article, this formulation can be slightly simplified. This section first details the topology of the considered class of networks, then the traditional MINLP formulation adapted to this class.

Notation and variables 195

In the considered class of water networks, illustrated by Figure 1, a pumping station is connected to a pipe network in which water is transported to the elevated water tanks. These can naturally be represented by a multi-arc K connected to an arborescence, i.e. a directed rooted 200 tree, (J, L):

• each arc k ∈ K models one pump and links the inlet reservoir node r to the outlet station node s;

• the set J of nodes in the arborescence includes the station node s (the root), the tanks J T (the leaf 205 nodes), and the pipe junctions J J (the internal nodes);

• the set L of directed arcs consists of the pipes L P and the flow control valves L V . In our context, there is exactly one valve before each tank. In the arborescence (J, L), the valves are then represented by 210 the terminal arcs to the leaves (from a junction to a tank), while the pipes are represented by the internal arcs (from the station node or between two junctions).

The planning horizon is one day divided into T = 24

215
hourly periods. The mathematical model involves two sets of decision variables: the activity indicator

x kt ∈ {0, 1} of pump k ∈ K (i.e. x kt = 1 if k is on at time t, x kt = 0 otherwise) and the flow rate q lt ≥ 0 (in m 3 /h) passing in a directed 220 link l ∈ L or pump l ∈ K in time period t ∈ [1, T ].
Since each time period t lasts 1 hour, q lt is also the volume (in m 3 ) of water flowing through element l during period t.

The model also involves one set of implied variables: the hydraulic head (or node potential) h jt ≥ 0 (in m) at node 225 j ∈ J in time period t ∈ [1, T ]. This is conventionally measured as the sum of the geographical elevation and the water pressure head at the node.

J nodes L arcs r reservoir node K pumps s ∈ J station node L P ⊂ L pipes J T ⊂ J tank nodes L V ⊂ L valves J J ⊂ J internal nodes t ∈ [1, T ] time periods x kt ∈ {0, 1} activity indicator of pump k ∈ K in period t q lt ≥ 0 flow/volume through l ∈ L ∪ K in period t h jt ≥ 0 hydraulic head at node j ∈ J in period t

Potential-flow coupling

The physical behavior of a pump or a pipe is accurately 230 modeled by a real function of the instantaneous flow rate passing through it:

Pipes. Darcy-Weisbach is the preferred law to describe the pressure loss along a pipe ij ∈ L P (as a function Φ ij of the flow q) in water networks [START_REF] Eck | Valve placement in water networks: Mixed-integer non-linear optimization with quadratic pipe friction[END_REF]. This equation relies on a coefficient, the friction factor, which is itself an implicit function of the flow [START_REF] Geißler | 845 Mixed integer linear models for the optimization of dynamical transport networks[END_REF]. To obtain a tractable expression within a mathematical model, this coefficient is generally approximated using the Prandtl-Kármán formula for fully turbulent flows (with a Reynolds number Re greater than 10 6 ) as in [START_REF] Verleye | Optimising production and distribution operations in large water supply networks: A piecewise linear optimisation approach[END_REF][START_REF] Gleixner | Towards glob-860 ally optimal operation of water supply networks, Numerical Algebra[END_REF] or using the Swamee-Jain formula for turbulent flows (Re > 4000) as in the modeling software EPANET [START_REF] Rossman | EPANET Users Manu[END_REF]. To keep our model both accurate and tractable, and to make it compatible with EPANET for checking the feasibility of the solutions, we follow Eck et al. [START_REF] Eck | Valve placement in water networks: Mixed-integer non-linear optimization with quadratic pipe friction[END_REF] and approximate the latter expression as a quadratic function

Φ ij (q) = A ij q + B ij q 2
over range 4000 ≤ Re ≤ 10 5 using a least squared regression method. For the two considered classes of pumps, the pressure discharge can accurately be approximated by a quadratic function (as in [START_REF] D'ambrosio | Mathematical programming techniques in water network optimization[END_REF]) of the flow:

Ψ k (q) = A k -B k q 2 .
Rather than expressing the power consumption as a function of pressure, flow and efficiency, we adopt a linear approximation: Γ k (x, q) = P 0 k x + P k q, where q is the flow rate through the pump and x the activity indicator of the pump. Note that, by definition 1 , the flow through a pump is null (q kt = 0) if and only if the pump is inactive (x kt = 0). The contribution to the cost is then null (Γ k (x kt , q kt ) = 0) when pump k is inactive at time t. Such a linear fit has been considered in previous 1 this definition is enforced by Constraints [START_REF] Ormsbee | Optimal control of water supply pumping systems[END_REF] in our model.

works (e.g., [START_REF] Geißler | 845 Mixed integer linear models for the optimization of dynamical transport networks[END_REF][START_REF] Verleye | Optimising production and distribution operations in large water supply networks: A piecewise linear optimisation approach[END_REF][START_REF] Skworcow | Pump schedules optimisation with pressure aspects in complex large-scale water distribution systems[END_REF]). In our case, it provides a very 245 good approximation as illustrated by the graphs on the right of Figure 2.

Finally, the pump efficiency is the ratio between the hydraulic power ρgqΨ k (q) (where ρ = 1000 kg/m 3 the water density and g = 9.81 m/s 2 the gravitational acceleration) 250 and the power consumed. Again, this quantity denoted η k (q) is well defined by a quadratic function (as in [START_REF] Verleye | Optimising production and distribution operations in large water supply networks: A piecewise linear optimisation approach[END_REF]).

Mathematical model

Given the notation summarized in Table 1, the Pump Scheduling Problem can be formulated as follows:

z = min t∈T k∈K C t Γ k (x kt , q kt ) (1) s.t. ij∈L∪K q ijt = ji∈L q jit , ∀t, j ∈ {s} ∪ J J (2) ij∈L q ijt = S j (h jt -h j(t-1) ) + D jt , ∀t, j ∈ J T (3)
h j0 = H 0 j , ∀j ∈ J T ( 4
)
H min j ≤ h jt ≤ H max j , ∀t, j ∈ J T (5) h it -h jt ≥ 0, ∀t, ij ∈ L V (6)
Q min k x kt ≤ q kt ≤ Q max k x kt , ∀t, k ∈ K (7) (h st -h rt -Ψ k (q kt ))x kt = 0, ∀t, k ∈ K (8) h it -h jt = Φ ij (q ijt ), ∀t, ij ∈ L P (9)
Objective (1) minimizes the total energy cost of pumping where C t denotes the electricity tariff (in e/kW h). Flow Constraints ( 6), ( 8) and ( 9) model the potential difference along, respectively, valves (as in [START_REF] Skworcow | Pump schedules optimisation with pressure aspects in complex large-scale water distribution systems[END_REF]), pumps and pipes. Finally, Constraints ( 7) and ( 8) model the activity of the pumps: Constraints (7) enforce the flow to be null when the pump is inactive, and within some positive interval 8) enforce a fixed potential difference between the station inlet r and outlet s only if at least one pump is active. If all pumps are off (x kt = 0, ∀k ∈ K), the hydraulic heads at nodes r and s are left unrelated.

270 [Q min k , Q max k ] otherwise. Constraints (

Towards a tractable model

Due to its non-convexities, model (1-9) above is likely hard to resolve at optimality. We show in this section how to exploit the specificities of the considered class of networks to turn this model into a tractable convex MINLP.

One-direction flow

Compared to the general formulation, the MINLP above presents a first simplification that is valid for any branched networks: pipes and valves allow for flows in only one direction. In the general case, either the flow variables q ijt on pipes are not restricted in sign (as in [START_REF] Gleixner | Towards glob-860 ally optimal operation of water supply networks, Numerical Algebra[END_REF]), leading to nonsmooth and nonconvex constraints [START_REF] Price | Iterative linearization scheme for convex nonlinear equations: Application to optimal operation of water distribution systems[END_REF], or the pipes are duplicated (q ijt ≥ 0, q jit ≥ 0) in the model and linked by additional binary variables (with x ijt + x jit ≤ 1) to model the flow direction (as in [START_REF] Eck | Valve placement in water networks: Mixed-integer non-linear optimization with quadratic pipe friction[END_REF]). For valves, equivalently, ei-290 ther a nonconvex constraint [START_REF] Burgschweiger | Optimization models for operative planning in drinking water networks[END_REF] or an additional binary variable [START_REF] Gleixner | Towards glob-860 ally optimal operation of water supply networks, Numerical Algebra[END_REF] is required in place of a linear constraint [START_REF] Van Staden | A model predictive control strategy for load shifting in a water pumping scheme with maximum demand charges[END_REF] to ensure that the valve pressure decrease is consistent with the direction of the flow.

Symmetric pumps

A pumping station commonly contains pumps with similar characteristics and ages. The set K of pumps can then be partitioned into classes C(K). Since all of the active pumps in a class work on the same operating point according to [START_REF] Nicklow | The ASCE task committee on evolutionary computation in environmental and water resources engineering. State of the art of genetic algorithms and beyond in water resources planning and management[END_REF], these pumps can be aggregated in the model. One way of doing this [START_REF] Skworcow | Pump schedules optimisation with pressure aspects in complex large-scale water distribution systems[END_REF] is to model the number of active pumps as an integer variable instead of the binary activity status x kt of each individual pump. This reduces the number of head-flow constraints [START_REF] Nicklow | The ASCE task committee on evolutionary computation in environmental and water resources engineering. State of the art of genetic algorithms and beyond in water resources planning and management[END_REF] to only one per pump class. An alternative [START_REF] Gleixner | Towards glob-860 ally optimal operation of water supply networks, Numerical Algebra[END_REF] is to order the pumps in a class arbitrarily and enforce the following symmetry breaking constraints:

x (k+1)t ≤ x kt , ∀t, C ∈ C(K), k, k + 1 ∈ C (10)
The benefits of aggregation or symmetry breaking are twofold. First, it helps the computation of the MINLP by reducing both its size and the size of its search space (from 2 |K| to C∈C(K) |C| + 1); note that in our experiments, the two approaches resulted in similar improvements. Second, it decomposes the operational problem into two steps: the MINLP determines how many pumps per class to activate, while the operators of the pumping station determine which physical pumps to select. Based on their practical knowledge of the network elements, they can more naturally and efficiently enforce operational constraints and preferences regarding pump lifetimes and maintenance needs. For example, they can balance the activity across the pumps in order to limit the number of switches and the inactivity period of each pump per day.

Convex potential-flow coupling constraints

As suggested by Gleixner et al. [START_REF] Gleixner | Towards glob-860 ally optimal operation of water supply networks, Numerical Algebra[END_REF], if a pipe ij ∈ L P and a pressure-reducing valve l ∈ L V are connected with zero demand at junction j ∈ J, then the pressure loss equation [START_REF] Price | Iterative linearization scheme for convex nonlinear equations: Application to optimal operation of water distribution systems[END_REF] for pipe ij can be relaxed into a lower bound:

h it -h jt ≥ Φ ij (q ijt ), ∀t, ij ∈ L P (9')
Indeed, the valve can dissipate any extra amount h it -Φ ij (q ijt ) -h jt ≥ 0 of pressure at junction j.

In our case, a valve precedes each water tower and no water demand occurs at junctions. This allows us to relax 315 the equality in all constraints (9) because on each path connecting the pumping station to a water tower, the flow ultimately goes through a valve. The valve then compensates for the artificial higher pressure losses along the path.

For the same reason, we propose to relax the potentialflow coupling equations [START_REF] Nicklow | The ASCE task committee on evolutionary computation in environmental and water resources engineering. State of the art of genetic algorithms and beyond in water resources planning and management[END_REF] for pumps as

(h st -h rt -Ψ k (q kt ))x kt ≤ 0, ∀t, k ∈ K (8')
to allow an artificial lower pressure h st at the station out-320 let s. Again, we assume that the valves will dissipate the actual higher pressures at the junctions, when the operating point of each pump k is pushed up on the head-flow curve Ψ k .

The following two propositions show that the relaxed 325 MINLP (P ) : (1 -7, 8 , 9 ) is equivalent to the original MINLP (P) : (1 -7, 8, 9) on the condition that the slopes of the power functions of the pumps are all equal. More precisely, and of practical interest, we show how to build a feasible solution for the original problem (P) given a re-330 laxed solution of (P ) by balancing the flow through the pumps (to push up the operating point of the pumps on their respective head-flow curves) and increasing the pressure progressively from the station to the dissipating valves (to propagate the actual lower pressure losses).

335

In Proposition 1, we consider the case in which all pumps are identical. We prove that both solutions cost the same: the derived solution is then optimal for (P) if the relaxed solution is optimal for (P ). Proposition 2 addresses the general case in which pumps of different classes can 340 be active at the same time. The cost of the two solutions may then differ by a factor of max k∈K P k -min k∈K P k in the worst case, where P k is the slope of the power function Γ k of each pump k. The conversion leads thus to the optimum of (P) if the relaxed solution is optimal for (P ) and 345 if all slopes P k are identical. Otherwise, the derived solution may not be optimal for (P) but the optimum of (P ) provides a lower bound for (P) and thus an estimation of the optimal gap of the derived solution.

Proposition 1. Assume all pumps in K are identical and 350 let denote Ψ k = Ψ and Γ k = Γ for all pumps k ∈ K. For each solution (x, q , h ) feasible for the relaxed program (P ), there exists a solution (x, q, h) feasible for the original program (P) and at same cost.

Proof: We build a solution (x, q, h) that coincides with 355 (x, q , h ) with the possible exception of the values of the hydraulic pressures h jt at the junctions j ∈ {s} ∪ J J and the values of the flows q kt through the active pumps k ∈ K.

First, the flow is evenly spread between the active pumps at each time step t ∈ T according to: q kt = q t for all k ∈ K 360 such that x kt = 1 with q t = k∈K q kt k∈K x kt . The pressure at the station outlet is set to h st = Ψ(q t ) (using r as the reference point, we assume h rt = 0 w.l.o.g.). Since Ψ is concave,

then h st ≥ k∈K Ψ(q kt ) k∈K x kt ≥ min k∈K|x k =1 Ψ(q kt ) ≥ h st .
each time step t ∈ T recursively along the paths from station s by: h jt = h it -Φ(q ijt ) with q ijt = q ijt for all ij ∈ L P . By induction from s, we prove that h jt ≥ h jt for all j ∈ {s} ∪ J J : assuming h it ≥ h it for some pipe ij ∈ L p , then h jt ≥ h it -Φ(q ijt ) ≥ h jt . As a consequence, (x, q, h) satisfies constraints ( 6) for all water towers j ∈ J T , valves ij ∈ L V and time steps t ∈ T since h it -h jt = h it -h jt ≥ h it -h jt ≥ 0. All other constraints of (P) are also satisfied by construction of (x, q, h).

Finally, since Γ is linear, both solutions cost the same:

z = t∈T k∈K C t Γ(q kt ) = t∈T k∈K C t Γ(q kt ) = z .
Proposition 2. For each solution (x, q , h ) feasible for the relaxed program (P ), there exists a solution (x, q, h) feasible for the original program (P) such that the cost difference is bounded from above by a factor of (max k∈K P kmin k∈K P k ).

Proof: The derived solution (x, q, h) is built as in Proposition 1 except for the distribution of the total flow through the pumps at each time step t.

385

The flow conveyed by the active pumps of a given class C ∈ C(K) can be evenly balanced as above. Let n Ct = k∈C x kt be the number of active pumps in class C and q Ct = k∈C q kt n Ct be the flow through each active pump of C if evenly balanced, then the total flow at time t is 390 C∈C(K) n Ct q Ct . The flow must now be balanced among the pump classes such that Ψ C (q Ct ) = Ψ C (q C t ) for all classes C, C ∈ C(K), and such that total flow is conserved, i.e. C∈C(K) n Ct (q Ctq Ct ) = 0. This system of |C(K)| equations and variables q Ct admits exactly one solution as the pressure discharge functions Ψ C are bijective on the half-plane q ≥ 0.

Let h st = Ψ C (q Ct ) for any c ∈ C(K). The proof of Proposition 1 applies to establish the feasibility of solution (x, q, h) if we prove that h st ≥ h st : Since C∈C(K) n Ct (q Ct -400 q Ct ) = 0, there is at least one class C ∈ C(K) such that q Ct ≤ q Ct . Given (8') and that Ψ C is both concave and decreasing, then

h st = Ψ C (q Ct ) ≥ Ψ C (q Ct ) ≥ h st .
Finally, the cost difference at time t between the two solutions is ∆C t = C∈C(K) C t n Ct P C (q Ct -q Ct ). Let the set of classes C ∈ C(K) be divided in C(K) + and C(K) - depending respectively on whether the flow is increased (q Ct ≥ q Ct ) or decreased (q Ct < q Ct ), and let P + and P -be respectively the highest and lowest slope of power functions among the classes C ∈ C(K). Then, ∆C t ≤ 410 C t P + C∈C(K) + n Ct (q Ct -q Ct )+C t P - C∈C(K) -n Ct (q Ctq Ct ) = C t (P + -P -) C∈C(K) + n Ct (q Ct -q Ct ). The total cost difference is thus bounded from above by a factor of (P + -P -).

Finally, the potential-flow coupling equations (8') can be reformulated as quadratic inequalities with the help of big-M values:

h st -h rt ≤ Ψ k (q kt ) + M (1 -x kt ), ∀t, k ∈ K (8")
By choosing a sufficiently large value for M (we use M = 415 max k∈K A k -min k∈K A k ), this constraint enforces h sth rt ≤ Ψ k (q kt ) when pump k is active (x kt = 1) and leaves h st and h rt unrelated otherwise.

The resulting model (P ) : (1 -7, 8 , 9 , 10) belongs to the class of convex Mixed Integer Quadratically Con-420 strained Programs (MIQCP), and thus has in practice a lower computational complexity than (P). Indeed, it can be handled by efficient solution algorithms capable of reaching global optimality, such as the sophisticated Branch-and-Bound and Outer Approximation algorithms 425 available in numerous general-purpose solvers. Furthermore, in order to derive a feasible solution for (P) with a performance guarantee, (P ) does not need to be solved to optimality, but only to provide at least one feasible solution and a lower bound. In the experiments we present 430 below, we run the solver for a limited amount of time, enough to compute this information.

Application to the case study

Study case and experimental set

In this section, we present the results of our experi-435 ments conducted on the medium-size network FRD depicted in Figure 1: 16 water towers with valves are connected by 31 pipes to a single pumping station equipped with 4 pumps of class KSB150 and 2 pumps of class KSB100 (the pump characteristics are depicted in Figure 2).

440

FRD is a real drinking water supply system located in a rural zone of France. It serves an area of about 350km 2 with a marked topography: the source is at elevation of +40m while the highest reservoir is at +131m. The system annually supplies about 3 million m 3 of water, with 445 twice as much in summer compared to winter. Currently the pumping station is operated by an in-house heuristic solution (this common strategy [START_REF] Feldman | Aspects of energy efficiency in water supply systems[END_REF] is described in section 5.3.1), with an annual consumption of about 1.3GW h. We have access to one year of historical data on the aggre-450 gated water demand of the end-consumers at each water tower. Data consist of measurements 10 or 15 minutes apart that we average out into a 1-hour step. The electricity tariff scheme we consider is the green tariff A5 Base Very Long Uses provided by the French utility company 455 EDF [START_REF]Arrêté du 30 octobre 2014 relatif aux tarifs de vente de l'électricité[END_REF]. This day/night pricing has two different tariffs for summer (April to October) and winter (November to March). Furthermore, a peak tariff applies during four hours a day (two in the morning and two in the evening) from December to February, and night pricing applies all 460 day Sunday.

We empirically evaluated our approach on 365 different instances, each corresponding to a day of the historical data. The water demands -and thus the instances -are highly contrasted depending on the season; in particular, 465 demand increases significantly during the summer tourist season. A 'day' starts and ends at 22h, which marks the beginning of the night pricing, and the water tower volumes are set to their minimum level at the beginning of the day. Indeed, to minimize the total energy cost in practice, the water towers are emptied at the end of the day and filled at the beginning of the next day, i.e. during the night pricing period.

The solution approach is as follows: as explained in Section 3, curves Γ C and Ψ C for the two classes of pumps are fitted from the series of operating points provided by the constructor and the pressure loss function Φ p is approximated using a quadratic function. Model (P ) is then built with the Gurobi Python API and solved by the default MIQCP solver of Gurobi v5.6.3 using 8 threads of a quad-core processor 8GB RAM Intel Core i7 CPU 3720 @ 2.66 GHz running Windows 7 (64 bits). Note that this solver is also freely available on the NEOS server [START_REF]The NEOS server[END_REF] through different interfaces. We allot 60 seconds to the computation time after which the solver returns the best lower bound and the best solution found. The solution is then converted into a feasible solution to the original problem (P) according to the construction of Proposition 2. Note that, with only two different classes of pumps, the redistribution of the flows q Ct through the pump classes can be easily computed by analytically solving the system of two quadratic equations.

Validating the approach

The experiments were conducted on the 365 instances, but to investigate the results more in depth, we exhibited two representative days: a winter day and a summer day. To validate our approach on our case study, we answered the following questions:

May the approach compute good solutions in short times ? For all 365 instances, Gurobi found a lower bound and at least one feasible solution for the convex relaxation (P ) in less than 1 minute. The process to convert the best relaxed solution computed in 1 minute into a feasible solution for (P) is immediate. We estimated for each day d = 1, . . . , 365 the gap G d between the cost of the converted solution and the lower bound of (P ) as a ratio of the lower bound. The mean gap 1 365 d G d and the maximal gap max d G d are respectively 3.0% and 8.2%. The gap includes the surplus cost arising from the conversion of the solution to (P ) into a feasible solution to (P). This conversion cost represents only 0.17% of the overall gap: this is less than e0.01 per day on average.

Is the convex relaxation required ?

We directly solved the non-convex formulation (P) with the symmetry breaking constraints [START_REF] Feldman | Aspects of energy efficiency in water supply systems[END_REF] under nine state-of-the art MINLP solvers available on the NEOS server [START_REF]The NEOS server[END_REF]. Three solvers (Dicopt, Lindo Global, SBB) were unable to process the model because of its size. Four solvers (Bonmin, Couenne, Scip, Knitro) were unable to compute any feasible solution in less than 2 hours. The two global solvers AlphaECP and Baron computed good solutions in reasonable time for some instances but not for all. Figure 3 illustrates the behavior of these solvers compared to our approach on the typical winter day (top) and summer day (bottom). On the easier (winter) instance, our approach found the best 525 solution in 2 seconds, at 2% from the best lower bound LB computed by our approach. AlphaECP was also able to compute instantly a solution, but of inferior quality at 9% from LB. The solution did not improve before 30 seconds and met our best solution after more than 1 hour.

530

Note that Baron and Scip computed good first solutions but, respectively in 10 minutes and 5 hours. On the hard (summer) instance, none of the non-convex solvers were able to compute at least one solution in less than 5 minutes. The best solutions found after 2 hours by AlphaECP 535 and Baron were respectively at 39% and 20% from LB. Our approach found its best solution at 7% from LB in less than 20 seconds. These experiments show that a direct non-convex solution approach is not practicable in this application case.

540

Is the approach robust to increasing demands ? We now detail our approach behavior on the two representative days. On the winter day d W , the overall water demand is low (6110 m 3 ). For this instance and 229 similar others, only the two smallest pumps are activated in the computed 545 solution. The solution for d W is depicted in Figure 4b. Its optimality gap is G d W = 2%. According to Proposition 1, the conversion cost from (P ) and (P) is zero.

On the typical summer day d S , the water demand is doubled (13505 m 3 ). For this instance and 134 others, the 550 largest pumps are also activated, supplying about 20% of the total flow. We observe that the solution optimality gap is greater when the daily flow is higher: for the day under study, the gap is G d S = 7%. The conversion cost is negligible (e0.002) for two reasons. On the one hand, 555 the slopes of the power functions of the small and large pumps are close. On the other hand, the largest pumps are active only on few time steps. For 51 of the 135 summer instances, the conversion cost is even zero because, at time steps where flows need to be rebalanced, the active pumps are all of the same class.

Are the computed solutions physically viable ? We tested the solutions with EPANET. To match the hypothesis of this simulator, we made two slight changes to our model. First, we enforced the pressure loss functions to overestimate the pipe losses, i.e. we computed A ij and B ij with the constraint that Ψ ij (q) is greater than the theoretical pipe losses for the range of possible flow values. Secondly, we replaced constraints (6) on valves with

h it -h jt ≥ Q M ax ij /S j , ∀t, ij ∈ L V (6') given Q M ax ij
an upper bound on inflows. This forces the water tower head to increase when inflows exceed demand. With these two modifications, we observed that our solutions are fully compatible with EPANET for both the winter day and the summer day.

Evaluating the gains achieved

In this section, we estimate the annual expected benefits arising from optimizing the pump schedule. To this end, we simulated the strategy currently implemented in the FRD network. Variants of this strategy are common in water utilities [START_REF] Feldman | Aspects of energy efficiency in water supply systems[END_REF].

Current strategy

In practice, the operators manage the network by opening and closing gate valves at the water towers. Additionally, flow control valves fix the amount of water entering each tower. This amount results from the trade-off between securing the water supply and restricting pressure losses; it is only modified twice a year to correspond to seasonal demand. At each valve status change, if the induced discharge pressure decreases or increases sufficiently, the operators select pumps to switch on/off. In our simulation, the time is discretized into 15-minute periods, and the combination of pumps to activate in a given period is selected so as to satisfy the required total flow and minimize power consumption.

The operators decide to open or close a valve at a given time when the current volume of the water tower reaches thresholds: the valve is open (resp. closed) if the volume is within or passes beyond the lower (resp. upper) threshold and remains open (resp. closed) until it reaches the upper (resp. lower) threshold. To compute these thresholds, the operators take into account the electricity tariff and, very approximately, the future demand: the range allowed between the thresholds has been estimated through practice to avoid over-frequent changes in valve status, and the threshold levels depend on the tariff period. Precisely, at night, the upper threshold is set to the physical capacity of the tower and the lower threshold is set slightly below this value in order to provoke pumping at the lowest tariff 600 and store for the day ahead. Conversely, during the day, the lower threshold is set at the minimum water volume required in case of emergency and the upper threshold is set slightly above this value in order to limit pumping at the highest tariff.

Estimated benefits

To estimate and explain the energy cost savings, we consider the alternative formula of the cost as a function of the flow, pressure and efficiency2 , denoted as z d for a day d = 1, . . . , 365 :

z d = t∈T k∈K q ktd ρgΨ k (q ktd ) η k (q ktd ) C td = t∈T C td k∈K q ktd Ψ k (q ktd ) η k (q ktd ) t∈T k∈K q ktd Ψ k (q ktd ) η k (q ktd ) C d t∈T k∈K q ktd Ψ k (q ktd ) η k (q ktd ) t∈T k∈K q ktd Ψ k (q ktd ) η -1 d × ρg t∈T k∈K q ktd Ψ k (q ktd ) t∈T k∈K q ktd Ψ d t∈T k∈K q ktd Q d
We then investigated the three factors of cost: the mean discharge pressure Ψ d needed in order to supply the water demand, the mean efficiency η d of the pumping station, and the mean cost C d of the electricity consumed. Thus, the significant 16.9% savings on the annual energy bill resulting from MINLP compared to BAU can be explained by at least one of these factors: better alignment with the power prices or lower power consumptionitself induced by lower discharge pressure or better pump 620 efficiency. Mean power cost. The alignment with the simple day/night tariff is comparable for the two strategies with an advantage for MINLP : the share of electricity consumed during the night is 17.4% greater for MINLP than for BAU and the share for the 4 hour-peak period is 96.6% lower. The mean cost of electricity decreases from e45.0/kWh for BAU to e43.1/kWh for MINLP . The price information is actually a key element of the scheduling decisions in both strategies, including in BAU where the fixed thresholds force the water towers to stay at their higher level during the night and at their lower level during the day.

Mean discharge pressure. MINLP provides a substantial reduction (13.3%) in annual energy consumption, driven by a 21.9% decrease in the mean discharge pressure. This lower pressure is clearly visible for the sample day depicted in Figures 4a and4b (see green diagrams): the discharge pressure in BAU is frequently well above 110m which is the maximum pressure generated in MINLP. The reason for this is twofold. First, MINLP fills the water towers rationally over the time horizon, while BAU enforces a rapid fill of the towers at some time steps, involving larger instantaneous flows and more pressure losses in the pipes.

In the illustrated example (see red diagrams), we observe that the filling of the water towers is well distributed over the off-peak period in MINLP , while in BAU , three large pumps send out a significant share of the water required during the first hour. Secondly, due to the binary status of the gate valves and the fixed amount of flow pass-ing through them, at any time BAU must select one set 650 of possible flow configurations from a finite number. Although the valves and pumps were jointly designed to fit the expected hydraulic conditions, any selected configuration of the valves may imply a substantial increase of the discharge pressure in order to reach the head-flow curves of 655 the pumps. This can be observed on Figure 4a at 8:15am: although the flow is unlikely to substantially increase the discharge pressure (the flow is less than 100 m 3 /h and serves a water tower located in a low-altitude region), the discharge pressure has to be raised to 140 meters in order 660 to intersect the head-flow curve of pump KSB100-1.

As an aside, the pressure decrease provided by MINLP also limits water leakages in the network. Indeed, water leakages are commonly modeled by an additional pressuredependent demand at nodes [START_REF] Eck | Valve placement in water networks: Mixed-integer non-linear optimization with quadratic pipe friction[END_REF][START_REF] Ulanicki | Open and closed loop pressure control for leakage reduction[END_REF] and the frequency 665 of pipe bursts also depends on the pressure [START_REF] Fantozzi | Including the effects of pressure management in calculations of short run economic leakage levels[END_REF]. Hence, minimizing the energy costs indirectly leads to positive externalities concerning leakages and maintenance costs.

Mean pump efficiency. The last factor we studied was the pump efficiency. We observed that the annual mean ef-670 ficiency was reduced from 74.2% in BAU to 67.2% in MINLP . Figure 5 illustrates this gap on the winter day, by comparing the occurrence of the pump operating points of BAU (blue) and MINLP (yellow): MINLP chooses to move away from the best efficiency flow in order to reduce 675 the discharge pressure, while the operating points in BAU are properly distributed all along the curves. 4a and4b, computed respectively with strategies BAU and MINLP : cumulated (bottom graphs), and projected on the head-flow curve (in red) and the efficiency-flow curve (in green).

Hence, MINLP significantly reduces the energy costs compared to BAU , despite operating the pumps at lower efficiency levels. It is worth noting that we could expect more savings if the pumping station were designed in compliance with the MINLP approach, e.g. by righttranslating the efficiency-flow curves of the pumps. To date, due to the limitation of the solvers, active network elements (like pumps, valves and tanks) are disregarded in the optimal design of water networks [START_REF] D'ambrosio | Mathematical programming techniques in water network optimization[END_REF]. However, our results call for a study of a systemic approach that takes into account both the investment and the operational costs.

Load Shifting

The previous experiments validated the MINLP ap-690 proach and illustrated the benefits in place of the current pumping strategy. Additionally, we estimated the robustness of MINLP to highly dynamic electricity tariffs.

While BAU is mainly driven by water volume thresholds, which take into account half-day-based pricing, MINLP is 695 theoretically compatible with hourly-based pricing. This type of strategy could thus be envisaged to implement load shifting at the pumping station, that is to take advantage of the water storage capacities in order to benefit from, and to contribute to, flexible power market conditions. To certify this usage, we experimented the same benchmark against the electricity price provided by the European Electricity Index [START_REF]ELIX european electricity[END_REF]. Figure 4c depicts the solution obtained by MINLP for a sample day. As expected, dynamic pricing does not hinder the capacity of the approach to compute near-optimal solutions in very short times: at least one solution was found in 60 seconds for all of the 365 instances, and, with respect to the day/night pricing, the value of the mean optimality gap is only slightly higher (3.6%) and the maximal gap is of the same order (9.6%). 710

Conclusion and Future Works

This paper proposed a new mathematical programming approach for pump scheduling in a common class of branched networks with one pumping station raising water to elevated tanks. Due to the presence of a flow control 715 valve at each water tower, the equality constraints arising both from the pipe pressure loss and the pump head-flow constraints can be relaxed, making the MINLP convex and thus tractable with general-purpose solvers. We then show how to easily derive a feasible solution to the original prob-720 lem and how to characterize its optimality gap.

We applied our approach to a real rural drinking water distribution network comprising 16 water towers, 31 pipes, and 6 pumps of 2 different classes. We experimented on 365 different instances, each corresponding to a day of his-725 torical water demand data. Using our approach, we were able to compute feasible solutions for all instances in less than 60 seconds with an average 3.0% optimality gap. The annual energy bill was 16.9% lower than a simulation of the current operational strategy based on water volume 730 thresholds. Finally, the approach succeeded in optimizing when subject to hourly-based electricity prices, showing that it is suitable for implementing demand-response, in contrast with the current strategy.

Our approach has two main advantages compared to 735 existing heuristic solutions: (1) it is able, in a very short time and for networks of a realistic size, to produce optimal or near-optimal solutions with performance guarantee, and (2) based on an adaptive model that can be handled by any appropriate black-box solver, it requires almost no 740 software development and is easily reproducible.

In future works, our first aim will be to generalize the approach to other water networks. As the convex model remains valid, we need to figure out ways to build feasible solutions from the relaxation and to characterize the 745 relaxation gap. Also, although a commercial solver can quickly solve the proposed convex program for the considered medium-size network, our approach still relies on solving an NP-hard problem. The limits of the approach in terms of network sizes and computation times should then 750 be investigated in order to address larger branched networks. Finally, the question is how to integrate different network elements (e.g. specific valves, serial pumping stations, demand at junction, variable speed pumps, loops) to the relaxation without altering its performance.

755

Another research focus is to expand our approach for operational control to the entire decision process driving the management of water networks at different time scales. For the short term, our approach can easily be extended to real-time scheduling and embedded in a closed-loop control 760 strategy [START_REF] Van Staden | A model predictive control strategy for load shifting in a water pumping scheme with maximum demand charges[END_REF]: as soon as the demand or the price diverges from the forecast, the model can be solved quickly, and the schedule can be repaired in line with the new solution. A practical question here is how the prediction errors deteriorate the solution quality. For the long term, our exper-765 imental results showed that significant benefits can be expected from tackling the whole optimization problem coupling design and scheduling in a systemic way; that is by simultaneously considering the operational costs with the investment and maintenance costs. MINLP approaches al-770 ready exist for both design and scheduling and are closely based on the same core modeling elements. One open question is how to integrate the dynamic elements (pumps and valves) at the scale of the network design problem.

  115

Figure 1 :

 1 Figure 1: FRD network : a pumping station operating 2 small and 4 large pumps in parallel and supplying 16 water towers equipped with flow control valves.

Pumps.

  The characteristics of a fixed speed pump k ∈ L are given by the manufacturer as a set of experimental points that we propose to interpolate as depicted in Figure 2.

Figure 2 :

 2 Figure 2: Pressure discharge (left, in m), efficiency (middle, in %), and power consumption (right, in kW ) of the two classes of pumps available at the FRD pumping station: experimental observations (in black) and fitted curves (in orange).

  255conservation is ensured by Constraints (2) at the station node s and at each junction j ∈ J J , and by Constraints (3) at each tank j ∈ J T . In Constraints (3), D jt denotes the demand (in m 3 ) in time period t ∈ [1, T ] and S j the crosssection (in m 2 ) of the cylindrical tank: as the water in the 260 tank is at atmospheric pressure[START_REF] Burgschweiger | Optimization models for operative planning in drinking water networks[END_REF], the hydraulic head h jt represents the elevation head of the filling level, i.e. the elevation of the tank plus the volume of water divided by S j . The volume of water in a tank, and thus the hydraulic head, are bounded at any time (5) and fixed at time 0 (4).

  265

  275

Figure 3 :

 3 Figure 3: Comparison of the solutions (quality and computation time) found by MINLP solvers (AlphaECP, Baron, Scip) and by our approach (Gurobi) on the winter (top) and summer (bottom) representative days. LB denotes the value of the best lower bound computed by our approach.

  605

Figure 4 :

 4 Figure 4: Comparison between the simulated strategy (BAU) and the optimized (MINLP) schedules on a winter day instance for two electricity tariffs (day/night and dynamic price). Top to bottom: in red, the schedule of the flow through the large (top four) and small (bottom two) pumps, the discharge pressure at node s, the electricity tariff.

Figure 5 :

 5 Figure 5: Operating points of the pumps for the two classes KSB150 and KSB100 in the daily schedules depicted in Figures4a and 4b, computed respectively with strategies BAU and MINLP : cumulated (bottom graphs), and projected on the head-flow curve (in red) and the efficiency-flow curve (in green).

Table 1 :

 1 Summary of notation

Table 2

 2

	610
	summarizes the values obtained by our approach (denoted
	MINLP ) compared to the simulated current strategy (de-
	noted BAU ) aggregated on the annual scale, while Fig-
	ures 4a and 4b illustrate them on the typical winter day.
		BAU	MINLP	∆M/B
	power cost power consumption d z d (ke) d Ψ d Q d η -1 d mean discharge pressure 1 365 d ψ d (m) (GWh) mean pump efficiency 1 365 d η d mean power cost 1 d C d (e/MWh) 365 night/total power consumption (%)	58.7 1.335 133.1 74.2% 45.0 46.1%	48.8 1.157 103.9 67.2% 43.1 54.1%	-16.9% -13.3% -21.9% -9.6% -4.2% +17.4%
	day/total power consumption (%)	51.0%	45.8%	-10.2%
	peak/total power consumption (%)	2.9%	0.1%	-96.6%

Table 2 :

 2 Annual energy cost and consumption: comparison between the simulated strategy (BAU ) and our approach (MINLP ). The deviation ∆M/B is computed as (MINLP -BAU )/BAU.

	615

tions[START_REF] Geißler | 845 Mixed integer linear models for the optimization of dynamical transport networks[END_REF][START_REF] Morsi | Mixed integer optimization of water supply networks[END_REF][START_REF] Verleye | Optimising production and distribution operations in large water supply networks: A piecewise linear optimisation approach[END_REF][START_REF] De La Perrière | Water planning and management: An extended model for the real-time pump scheduling problem[END_REF]; or the time-coupling constraints are

Assessed over the

instances, the cost difference between this formulation and the linear expression used as the objective (1) of the MINLP is around 3%.

Acknowledgments

This work was partially funded by the Artemis Arrowhead project under grant agreement number 332987 (Artemis/Ecsel Joint Undertaking, supported by the European Commission and French Public Authorities).