Julian Fey
email: julian.fey@fmf.uni-freiburg.de

Simon Procz

Michael K Schütz

Vincent Schoepff
email: vincent.schoepff@cea.fr

Frédérick Carrel

Juan S Useche

Alex Fauler

Michael Fiederle

CdTe sensor configurations for Robot assisted photon counting gamma camera

Keywords: CdTe, Photon counting, Timepix, Gamma camera, Spectroscopy

In this work, different Cadmium Telluride (CdTe) sensor configurations are assessed for the usage in a robot assisted portable gamma camera. In the first part, four CdTe sensors, with thickness of 0.45 mm, 1 mm, 2 mm and 3 mm and pixel sizes of 55 μm and 110 μm, are investigated regarding their spectroscopic performance. The photon counting detector Timepix1 is hereby used. The 3 mm CdTe sensor shows increase in count rate up to a factor of 1.25 compared to a 2 mm CdTe sensor, 1.84 compared to a 1 mm CdTe sensor and up to 2.71 compared to a 0.45 mm CdTe sensor in the case of 137 Cs. In the second part, the 3 mm CdTe sensor was implemented in a commercially available gamma camera, the iPIX. The system was integrated in the bomb disposal robot and tested in different scenarios. The integrated 3 mm CdTe detector measured 21.5 counts per second emitting from a 60 Co source with an activity of 2.8 ± 0.07 Gbq in 20-meters distance in an open environment. The acquisition time was 116 seconds. The angular resolution was sufficient for the user to localize the radioactive isotope inside the test structure.

Introduction

The threat to public safety posed by improvised explosive devices (IEDs) has increased significantly in recent years [START_REF] Van Ballegooij | The fight against terrorism[END_REF]. In the context of IED, one of the most important aspects, both in terms of real danger and media presence, is the threat of radiological dispersal devices (RDDs). The primary goal of RDDs is neither the damaging of infrastructure nor the injuring of people by the explosive effect, but rather, the contamination of the largest possible area by radioactive substances. The current development state of semiconductor sensors for the detection of both, hard X-ray and gamma-radiation allows to produce highly sensitive, compact gamma cameras. This work addresses the needs of the first responder services to this threat by researching and improving contactless, robot assisted methods for the investigation procedures. The aim is to provide end users with reliable, fast and safe means to investigate IEDs and evaluate the imposed threat.

The material as well as the thickness of the sensor are deciding factors of the performance of a gamma camera. High Z semiconductor materials, such as Cadmium Telluride (CdTe), Cadmium Zinc Telluride (CdZnTe) and Thallium bromide (TIBr), offer fast analyzation of the investigated isotope due to their high absorption efficiency. The absorption efficiency can be further improved by increasing the thickness of the semiconductor sensor. However, the effect of charge sharing is increased accordingly [2 -4]. Charge sharing describes the spreading of information from charge carriers over multiple pixels.

The effect is heavily depended on the photon energy of the event and the thickness of the material due to the drift of charge carriers through the sensor layer. The main driving forces are repulsion between charge carriers and diffusion. Thicker sensors archive a higher attenuation efficiency which leads to superior energy resolution to identify the isotope, especially in the high energy region. A smaller pixel size yields a better spatial resolution for localization of the isotope in the region of interest. A crucial aspect in the conception of a Gamma-or Compton camera lies in the dependence of absorption efficiency and the negative effect of the charge sharing on the spatial resolution.

In this work, four different ohmic sensor configurations were evaluated for the usage in a gamma camera. The material was supplied by Acrorad (Japan) while the flip-chip bonding was performed by the Freiburg Materials Research Centre (Germany). The 2 mm CdTe sensor was bonded to a Timepix1 readout chip with an increased pixel pitch of 110 μm. An increased pixel pitch (pp) counters the charge sharing effect on a hardware level at the detector side. The 0.45 mm, 1 mm and 3 mm CdTe detectors were connected to a Timepix1 readout with 55 μm pp. In this case, the 3 mm CdTe Timepix1 shows the highest effect of charge sharing but also provides the highest attenuation efficiency in photon energies above 100 keV (Fig. 1). Cluster analysis was performed for all detectors since all detector assemblies exhibit charge sharing regardless of configuration.

Materials and Methods

Cadmium Telluride

CdTe is an excellent semiconductor material for hard X-ray and gamma-ray detection. The high atomic number ( ≈ 50) and density (ρ = 5.854 g/cm³ at RT) provide superior attenuation efficiency above 30 keV compared to Si ( = 14) and GaAs ( ≈ 32) and the wide bandgap energy (Egap= 1.44 eV) permits the usage of the detector at room temperature without active cooling. The absorption efficiency can be enhanced furthermore by increasing the thickness of the material [START_REF] Fey | Investigations on performance and spectroscopic capabilities of a 3 mm CdTe Timepix detector[END_REF].

CdTe was considered and used as a semiconductor sensor for hybrid photon counting detectors (HPD) of the Medipix chipset family since the beginning of the Medipix collaboration. Previous work investigating the spectroscopic performance of a 1 mm CdTe sensor bump bonded to a Medipix2 chipset was presented by Maneuski et al. [START_REF] Maneuski | Imaging and spectroscopic performance studies of pixellated CdTe Timepix detector[END_REF], Ruat and Ponchut [START_REF] Ruat | Characterization of a Pixelated CdTe X-Ray Detector Using the Timepix Photon-Counting Readout Chip[END_REF] and Greiffenberg et al. [START_REF] Greiffenberg | Energy resolution and transport properties of CdTe-Timepix-Assemblies[END_REF].

Measurements regarding the characterization and spectroscopic performance of a Timepix1 chipset in combination with 2 mm CdTe [9, 10] and 3 mm CdTe [START_REF] Fey | Investigations on performance and spectroscopic capabilities of a 3 mm CdTe Timepix detector[END_REF] were shown recently. Cadmium Zinc Telluride (CdZnTe) is another promising semiconductor material in addition to CdTe. Recent work with CdZnTe sensors bump bonded to a Timepix1 chipset [START_REF] Tsigaridas | X-ray imaging with high-Z sensors for the ESRF-EBS Upgrade[END_REF][START_REF] Tsigaridas | Characterization of pixelated CdZnTe sensors using MAXIPIX[END_REF] and performance investigations [START_REF] Veale | Characterization of the Uniformity of High-Flux CdZnTe Material[END_REF][START_REF] Thomas | Characterisation of Redlen high-flux CdZnTe[END_REF] show encouraging results and are worth investigating in future work.

Timepix

The Timepix1 is a hybrid pixel detector (HPD) developed by the Medipix2 collaboration CERN [START_REF] Cudié | Design and characterization of 64K pixels chips working in single photon processing mode[END_REF]. It is fabricated in a CMOS 250 nm process, consists of 256 x 256 square pixels with a pixel pitch of 55 μm, resulting in a sensitive area of 1.98 cm 2 .The detector provides three different counting modes [START_REF] Cudié | Design and characterization of 64K pixels chips working in single photon processing mode[END_REF]: single photon counting mode (Medipix Mode), Time-Of-Arrival Mode (ToA) and Time-over-Threshold Mode (ToT). For this work, all measurements were carried out using the Time-over-Threshold Mode (ToT).

The incoming signal is, individually for each of the 65536 pixels, amplified by a charge sensitive preamplifier and compared to a 4-bit threshold in the discriminator circuit which is connected to a synchronization logic. If the signal exceeds the threshold, a 14-bit counter is incrementing the number of clock units as long as the signal is above the threshold and below 11810 counts in which case the overflow logic stops. The number of incrementing clock pulse is proportional to the input charge and therefore to the photon energy of the event due to a linear discharge of the preamplifier controlled by the IKrum current. This contrasts with single photon counting pixel detectors which discharge using a resistor resulting in exponential decrease of the signal hub [START_REF] Johnson | a single-photon counting x-ray detector[END_REF]. The energy threshold enables compensation of process variations.

Different semiconductor sensor materials like Si, CdTe, CdZnTe or GaAs can be connected to the Timepix ASIC and handled by the Krummenacher preamplifier which can process positive and negative charge inputs. The IKrum digital-to-analog-converter (DAC) (Figure 2) and has two functionalities: the compensation of leakage current and the discharge of the preamplifier feedback capacitor, which influences the duration of the return of the peak signal amplitude to the baseline. A

higher IKrum current leads to a faster discharge and thus shorter pulses which results in a reduction of pulse pileup but also to a reduction of signal amplitude, decreasing the energy resolution. The simulated dependence of the preamplifier output signal and width on the IKrum DAC settings are

shown in [START_REF] Hamann | Investigation of GaAs:Cr Timepix assemblies under high flux irradiation[END_REF]. in the case of electron collection IKrum/2, and the return of the signal to the baseline (adopted from [START_REF] Llopart | TimePix manual v1.0[END_REF][START_REF] Hamann | Investigation of GaAs:Cr Timepix assemblies under high flux irradiation[END_REF]).

The Timepix1 and its successor, the Timepix3, were initially aimed for particle detection.

However, in recent years these detectors are used in medical imaging [START_REF] Procz | X-ray and gamma imaging with Medipix and Timepix detectors in medical research[END_REF], space radiation dosimetry as well as in commercially available gamma cameras [START_REF] Amgarou | A comprehensive experimental characterization of the iPIX gamma imager[END_REF] and experimental Compton cameras [START_REF] Turecek | Application of Timepix3 based CdTe spectral sensitive photon counting detector for PET imaging[END_REF][START_REF] Amoyal | Development of a hybrid gamma imager for nuclear industry applications[END_REF] with 1 mm and 2 mm CdTe sensor configurations.

From GAMPIX prototype to iPIX gamma camera

In order to address the need of instrumentation allowing the detection and location of radioactive materials, CEA List has been developing gamma camera systems for over thirty years [START_REF] Gal | CARTOGAM -a portable gamma camera for remote localization of radioactive sources in nuclear facilities[END_REF]. These systems allow remotely visualizing gamma rays originating from radioactive emitters by superimposing a gamma image onto a visible image of the scene. Since the late 2000's, CEA List has been working on a new generation of gamma camera, which resulting prototype is called GAMPIX [START_REF] Gmar | GAMPIX: a new generation of gamma camera[END_REF]. This prototype relayed on the use of a Timepix1 detector equipped with a 1 mm thick CdTe sensor, coupled with a coded mask ensuring the spatial modulation of incoming gamma flux Figure 3 (a).

Coded aperture is a technique of great interest, which allows achieving a good angular resolution for gamma source reconstruction, while maintaining a high sensitivity of the system. GAMPIX gamma camera is based on MURA coded masks, consisting in a mathematical arrangement of holes, defining the rank of the coded mask, into a denal collimator [START_REF] Gottesman | New family of Binary Arrays for Coded Aperture Imaging[END_REF]. In the presence of radioactive sources, gamma-ray will go through holes and be stopped into denal, resulting in the projection of the MURA pattern on the gamma-sensitive sensor. The high pixelization of Timepix1-based detector allows accurately determining the location of gamma rays into the semi-conductor sensor, which ensures angular resolution of the overall imaging system. Indeed, the preferred pixel pitch of 55 µm ensures an achieved precision on pattern projection estimated around 0.20°, compared to 0.40° for a pixel pitch of 110 µm. This projection accuracy is a determining factor for the source location reconstruction. Moreover, the symmetry in the MURA pattern allows, by a 90° rotation of the coded mask, suppressing background noise, without any additional shielding. This feature, called anti-mask procedure, is of great interest in the case of measurements in areas subject to high gamma background while preserving lightness of the camera.

The GAMPIX prototype was validated in laboratory [START_REF] Amoyal | Metrological characterization of the GAMPIX gamma camera[END_REF] and deployed in various real environments in the field of nuclear industry and Homeland Security applications [START_REF] Carrel | GAMPIX: A new gamma imaging system for radiological safety and Homeland Security Purposes[END_REF]. Performance assessment of the gamma camera prototype showed that the system was able to perform imaging measurements in a field of view of 50° with an angular resolution ranging from 1.35° to 3.81° depending on the pattern of coded mask used. As presented in Table 1, the sensitivity of the prototype makes it possible to visualize sources in less than a minute over the whole range of energies encountered in the field of nuclear industry (i.e.

from 241 Am at 59.5 keV to 60 Co at 1.33 MeV) [START_REF] Carrel | GAMPIX: A new gamma imaging system for radiological safety and Homeland Security Purposes[END_REF].

Table 1. Sensitivity of GAMPIX gamma camera, according to the rank of the mask (R7 is for rank 7, R13

for rank 13, e2, e4, e8 indicates the corresponding thickness in mm) and to the nature of the radioactive source. The dose rate was measured in the vicinity of the gamma camera.

Radionuclide

Results

A sufficient bias was applied to all investigated detector assemblies during the presented measurements to ensure full charge collection efficiency. This was verified by measuring the count rate of the detector versus increasing sensor bias voltage under identical X-ray tube settings. The maximum charge collection efficiency is achieved when the count rate reaches saturation. The leakage current was measured before the experiments with no radiation present (Table 2). The high leakage current of the 0.45 mm CdTe is due to the effect of the processing of the thin material. A new dicing process was implemented which permits significant cleaner edges of the bonded sensor. The improved performance of this process can be seen in the low leakage currents of the 1 mm and 3 mm thick detector assemblies. In the case of the 2 mm sensor, the observed increase in leakage current is due to an earlier bonding date of the detector assembly. There is no surface preservation on the as-diced chip edge. Over time, degradation effects occur which lead to greater surface leakage currents. The DAC for leakage current compensation (IKrum) and the clock frequency were set to the same value for all assemblies to ensure comparable results regarding the spectroscopic performance (Table 2). All measurements were taken with a clock frequency of 9.6 MHz and IKrum 5 which corresponds to an IKrum current of roughly ~7.5 nA [START_REF] Cudié | Design and characterization of 64K pixels chips working in single photon processing mode[END_REF]. At IKrum 5, the mean pulse duration of a 100 keV photon is ~6 μs [START_REF] Llopart | TimePix manual v1.0[END_REF]. All detectors passed stability measurements with radiation present and continues operation in the range from 12 to 24 hours with the configurations shown in Table 2.

Even while performing cluster analysis, the charge sharing effect has a strong influence on the energy resolution of the sensor. This is caused by the limited rise time of the preamplifier in the pixel cell and the higher number of pixels marking the edge of the cluster, which have deposited energy of the event, but the threshold was not exceed. This, so called dead energy, has a greater impact on the energy resolution with increasing cluster size of the event since more surrounding pixels exhibit the dead energy effect.

Spectroscopic Performance

The common radioactive sources 137 Cesium ( 137 Cs), 152 Europium ( 152 Eu) and 60 Cobalt ( 60 Co) are investigated by the different detector assemblies in a laboratory environment. All detectors were calibrated by the noise edge method [START_REF] Procz | Optimization of Medipix-2 Threshold Masks for Spectroscopic X-Ray Imaging[END_REF] in Pixelman 2.2.3 [START_REF] Turecek | Pixelman: a multi-platform data acquisition and processing software package for Medipix2, Timepix and Medipix3 detectors[END_REF]. The acquisition time was determined for each setup individually to minimize event pile up. Since the thickest sensor exhibits the highest attenuation efficiency, a time was chosen with, on average, 5 events per frame for the 3 mm CdTe detector. All measurements were performed in a fixed timeframe of 1000 s and a fixed distance from source to detector of 60 mm. Cluster analysis was performed for each measurement which identifies counting pixels exhibiting charge sharing as well as separated fluorescence photons and contributes them to their original event. Gaussian fitting is used for the primary gamma ray line of each spectrum and the sigma (σ) is calculated with σ = w/2. Fig. 4, Fig. 5 and Fig. 6 show the measured spectra for 137 Cs, 152 Eu and 60 Co. The raw spectra were normed in ToT-value on the x-axis a distinguished feature in the spectrum were adjusted to the correlated energy because the TOT values differ per detector. The spectra are superimposed on each other to visualize the difference in performance. All spectra were smoothed by moving average of 10 and the count rate was normalized.

The relative absorption efficiency of all sensors was evaluated in comparison to the 3 mm CdTe sensor which exhibits the greatest attenuation efficiency. Previous work investigated the relative absorption efficiency of 1 mm and 0.45 mm CdTe sensors [START_REF] Fey | Investigations on performance and spectroscopic capabilities of a 3 mm CdTe Timepix detector[END_REF]. The data was corrected for the different number of inactive pixels per sensor caused by crystal inhomogeneities and degenerating effects like corner peel. The measurement error was calculated by Gaussian error propagation law.

3.1.1. Cesium-137 137 Cs is an isotope used in the medical environment for radiation therapy. Therefore, in the context of the evaluation of the performance of a gamma camera. The acquisition time was set to 10 ms for all detectors. The lab source has an activity of 3.58 ± 0.09 MBq at the time of the experiment. Figure 4 presents the measured spectra of 137 Cs with the different detector assemblies. A relative absorption efficiency of (0.370 ± 0.003) for 0.45 mm CdTe, (0.544 ± 0.002) for 1 mm CdTe and (0.801 ± 0.002) for 2 mm CdTe in comparison to the 3 mm CdTe was measured.

Europium-152

152 Eu is a gamma and beta emitter used in hard-gamma brachytherapy. At the time of the measurement, the lab source had an activity of 0.28 ± 0.01 MBq. The set acquisition time for all sensors was 32 ms, which was chosen to minimalize event pile up. Figure 5 shows the 152 Eu spectra measured by the four investigating detector configurations. For 152 Eu, the difference in absorption efficiency between the investigated sensors is less evident than for 137 Cs and 60 Co. This is due to the Sm L-lines at 40 keV (Sm Kα1) and 39 keV (Sm Kα2), originating from the 152 Eu decay chain, which have an absorption efficiency above 99 % for all investigated sensors.

In addition, the emitted beta particles of 152 Eu are absorbed in the first few μm of the sensor layer regardless of the thickness. However, multiple gamma ray lines at energies >200 keV contribute to the observed ratio in relative absorption efficiency compared to the 3 mm sensor of (0.637 ± 0.003) for 0.45 mm CdTe, (0.777 ± 0.003) for 1 mm CdTe and (0.944 ± 0.003) for 2 mm CdTe.

Cobalt-60

60 Co is used in medical radiotherapy as a radiation source and therefore common in most hospitals and a potential isotope in an IED. In this investigation, it is used to evaluate the spectroscopic performance of the detectors for high energy photons since 60 Co exhibits two characteristic gamma ray lines at 1173 keV and 1333 keV. At the time of the experiment, the used lab source had an activity of 6.55 ± 0.16 Mbq. The acquisition time was set to 5 ms. Figure 6 presents the measured spectra of 60 Co investigated by the different detector assemblies in this work. A significant reduction in overall count rate can be overserved for the 1 mm and 0.45 mm sensors with a relative absorption efficiency of with (0.238 ± 0.002) for 0.45 mm CdTe and (0.396 ± 0.004) for 1 mm CdTe compared to the 3 mm CdTe. The 2 mm CdTe sensor exhibits the same effect as the thinner sensors but to a lesser degree with a relative absorption efficiency compared to the 3 mm sensor of (0.641 ± 0.005).

Implementation and field testing

After investigating the spectroscopic performance in chapter 3.1 and 3.2, a field test campaign was initiated to evaluate the influence of a thicker sensor on the localization capabilities of the gamma camera. The 3 mm CdTe sensor has the highest attenuation efficiency of the investigated sensors and the 55 μm pixel pitch permits the necessary spatial and angular resolution for the localization of the threat. Therefore, a 3 mm CdTe Timepix1 detector assembly with 55 μm pixel pitch was implemented in the iPIX gamma camera. Figure 7 shows the camera mounted on the robot Telemaxx manufactured and provided by Telerob GmbH (Germany). A specific mounting solution of the gamma camera was provided by ZIPPERMAST GmbH (Germany) which mounts the iPIX on the robotic arm. Following this approach, it was possible to remain the manipulation capabilities of the robot on investigating objects while having a 360° radius for the FOV of the gamma camera.

NBC testing side

During the project, it was possible to carry out a field test of the modified iPIX at the NBC testing side of the Austrian army at the DABSCH Kaserne (Austria). A room was prepared for investigation with a simili, containing a radioactive source, placed inside. The complete system could be tested under realistic conditions outside the lab environment [START_REF] Czetina | Robot assisted analysis of suspicious objects in public spaces using CBRN sensors in combination with high-resolution LIDAR[END_REF]. Figure 8 (a) shows the superimposed image of the gamma camera with localized radiation source in a suitcase. 137 Cs is the investigated isotope with an activity of 36 ± 0.9 Gbq. The distance from source to detector was 1.7 meters with 218 recorded counts per second (Cps). Gbq is placed inside and localized by the 3 mm CdTe sensor with 218 Cps. (b) The spectrum which was recorded at the same position in 100s. The mask R7 -e4 (Table 1) was used for the measurement.

In comparison to the normalized 137 Cs spectrum measured in the lab, shown in Figure 4, the raw

ToT-values are presented on the x-axis in Figure 8 (b). A total measurement time of 100s was chosen since it represents a realistic time frame to investigate a suspicious object in case of an emergency situation regarding the potential isotope present. The overall spectrum is noisier which is a result of the lower event count. However, the gamma ray line at 662 keV can be distinguished.

Nuclear Powerplant test facility

A second field test was performed with the iPIX configuration described in section 3.3 at a nuclear power plant test site. A 60 Co source was placed inside a test structure which is 150 cm wide and 100 cm in height. Three tubes are in equal distance and the isotope is placed randomly in one arm. The source has an activity of 2.8 ± 0.07 Gbq which corresponds to a dose rate of ~1 mSv/h in 1 m distance. Table 3 shows the results of the measurements. Figure 9 presents the distance 25 meters (a) and 20 meters (b). A spectrum was recorded at 15 meters distance between the 60 Co source and detector shown in Figure 11. The total measurement time was 50 s due to on-site time restrictions. The spectrum is noisier than the spectrum in the lab environment. This is expected due to the distance of 15 meters to the source resulting in a low count rate. The information would not be feasible to identify the investigated isotope, however, the presence of high energy photons can be determined since the cluster analysis detects and sufficiently excludes alphas and muons from the evaluation [START_REF] Fey | Investigations on performance and spectroscopic capabilities of a 3 mm CdTe Timepix detector[END_REF].

Figure 11. 60 Co spectrum recorded during the measurement series at the nuclear power plant test facility. The distance between source and detector is 15 meters with a total measurement time of 50 s. The mask R7 -e4 (Table 1) was used for the measurement.

A test parkour was prepared inside the power plant which is designed to evaluate the maneuverability and manipulation capabilities of the robot as well as localization of hidden radioactive sources. Video 1 shows a short clip of the iPIX mounted on the Telemaxx at the parkour during ENRICH 2019 (Austria).

Supplementary Material Video

Video 1. This video shows the Telemaxx bomb disposal robot with the iPIX mounted on the side arm during ENRICH 2019 (Austria). The iPIX used the 3 mm CdTe detector configuration. The parkour is designed to test the maneuverability and manipulation possibilities of the robot. In addition, radioactive sources must be localized.

Comparison to related work in the field

The presented approach uses the small pixel pitch of the Timepix ASIC and cluster analysis to guarantee a high angular resolution of the gamma camera by means of the MURA coded aperture as well as to provide spectroscopic information. A comparable gamma camera design approach is the Caliste-HD [START_REF] Maier | Second generation of portable gamma camera based on Caliste CdTe hybrid technology[END_REF], which uses a CdTe hybrid detector with a 16 x 16 pixel matrix with 580 μm pixel pitch and MURA coded mask for the localization capabilities. The pixel size used in Caliste-HD is much higher than in the IPix to reduce charge sharing effects without performing a cluster analysis. The relation of sensor thickness to pixel pitch enables the evaluation of single pixel events only for spectroscopic investigation. This results in a superior energy resolution. The gamma ray line of 152 Eu at 121 keV is resolved with a sigma of 0.9 keV in comparisons to the sigma of 11 keV of the 3 mm CdTe Timepix used in this work. However, as both coded-aperture imaging systems use a MURA mask, the highly pixelated Timepix has an angular resolution ranging from 1.3° to 3.8°, depending on the isotope, while the Caliste-HD system is reporting an angular resolution of only 7° in the case of 152 Eu [START_REF] Maier | Second generation of portable gamma camera based on Caliste CdTe hybrid technology[END_REF].

Scintillator based gamma cameras show a 6.0° angular resolution at 356 keV with a LaCl3(Ce) scintillator consisting of 22 x 22 voxels using URA coded aperture [START_REF] Lee | A Dual Modality Gamma Camera Using LaCl3(Ce) Scintillator[END_REF] and 5.6° angular resolution with a CsI (Na) detector module consisting of 20 x 20 voxels using a MURA coded aperture [START_REF] Lee | Investigation of nuclear material using a compact modified uniformly redundant array gamma camera[END_REF].

The CsI (Na) compact hybrid gamma camera was used to investigate the position of UO2 pellets with 4.1 % enriched 235 U and an activity of 1.14 μCi [START_REF] Lee | Investigation of nuclear material using a compact modified uniformly redundant array gamma camera[END_REF]. The results were used to predict the minimum required detection time and count rate for monitoring different amounts of UO2 pellets at 10-meters distance between source and detector. A localization time of 199.4 s and a count rate of 30.9 cps was found feasible to use the CsI (Na) compact hybrid gamma camera in the proposed scenario. In this work, it was possible to localize a 60 Co source in 20-meters distance in 116 s with a significant lower count rate per second of 21.5 Cps. The estimated required time of localization in a time crucial situation agrees with the localization times stated in this work.

A single layer Compton camera is another camera concept which uses Compton scattering for isotope localization by reconstructing the scattering angle [START_REF] Turecek | Single layer Compton camera based on Timepix3 technology[END_REF]. The advantage to renounce the coded-aperture system enables the construction of smaller and lighter camera systems. However, since the localization principle relies on Compton scattering, the detection efficiency of the camera is declining with decreasing energy of the investigated isotope. Furthermore, the distance between source and camera has to be known for position identification.

Recently, a hybrid Compton camera was proposed, combing a classical Compton camera and a pinhole camera in a single detector system. This hybrid camera achieves a sigma of 3.5 keV at the gamma ray line of 137 Cs at 662 keV, in comparison to sigma 55 keV of the 3 mm Cdte Timepix, and an angular resolution of 10° [START_REF] Omata | Performance demonstration of a hybrid Compton camera with an active pinhole for wide-band X-ray and gamma-ray imaging[END_REF].

. Discussion and Conclusions

Four different CdTe sensor configurations were bonded on Timepix1 readout ASIC to evaluate the spectroscopic performance for the integration in a gamma camera used on a robotic platform for counter measurements for radioactive threats.

Energy resolution is systematically better with greater pixel size and thinner sensors due to less charge sharing present. The 2 mm CdTe 110 μm pp detector exhibits a sigma 34 keV for the 137 Cs gamma ray line at 662 keV compared to sigma 55 keV for the 3 mm CdTe 55 μm pp detector. A negative influencing factor on the energy resolution is the deposition of not registered energy around a cluster due to the threshold of the individual pixel cell. This effect increases with overall cluster size due to a greater number of surrounding pixels. However, the attenuation efficiency of a sensor is an important aspect for a gamma camera since it allows efficient charge collection for a fast localization of the radioactive source especially for high energy gamma rays which pose an extensive threat to public safety. For the common isotope 137 Cs, the 3 mm CdTe sensor shows increase in count rate up to a factor The relative absorption efficiency measured in section 3.2 allows the calculation of the estimate acquisition time of a 1 mm CdTe sensor in the same scenario. In 20-meters distance from the source, the 1 mm CdTe sensor would have a predicted 8.5 cps resulting in 293 s acquisition time. For the intended use of the gamma camera to assisted in situations where the public safety is concerned, the time difference in localization of the source is significant for immediate decisions on appropriate countermeasures. The angular resolution was sufficient for the user to localize the radioactive isotope inside the test structure without knowing the placement for this measurement. At a 25-meters distance between detector and source, it was not possible for the user to distinguish the placement of the source inside the test structure.

A teleoperated robot platform with manipulation capabilities in combination with a gamma camera can assist police, fire workers and personal from nuclear powerplants in hazard and security operations. It poses the prospect to position the user and personal in a safe distance from a radioactive threat while, at the same time, remain the possibilities of interaction with the investigated object.

Furthermore, in the frequent case of non-attacks, in which suspicious objects turn out to be harmless, the 3 mm CdTe sensor improves investigation speed significantly compared to the 1 mm CdTe sensor.

For future work, it is worth investigating even thicker CdTe sensors in combination with highly pixelated ASICs such as the Timepix1 as well as other high Z sensor materials such as CdZnTe and TiBr which have similar or superior attenuation efficiency as CdTe while are available in greater sensor thickness.

Figure 2 .

 2 Figure 2. Schematic of the Krummenacher circuit. The IKrum current controls the leakage current compensation,

Figure 3 .

 3 Figure 3. (a) Picture of GAMPIX second-generation gamma camera developed by CEA List; (b) Picture of iPIX gamma camera, industrialized version of GAMPIX and commercialized by MIRION Technologies.

Figure 4 .

 4 Figure 4. Spectra of 137 Cs investigated by four different detector configurations in this work. The sigma of the gamma ray line at 667 keV is shown.

Figure 5 .

 5 Figure 5. Spectra of 152 Eu investigated by four different detector configurations in this work. The sigma of the gamma ray line at 122 keV is shown.

Figure 6 .

 6 Figure 6. Spectra of 60 Co investigated by four different detector configurations in this work.

Figure 7 .

 7 Figure 7. Telemaxx by Telerob with the iPIX mounted on the arm. The mounting method does not negatively influence the manipulation ability of the robot while permitting at the same time a 360° radius for the position of the gamma camera.

Figure 8 (Figure 8 .

 88 Figure 8. (a) Investigation of a prepared smili at the NBC testing side. A 137 Cs source with an activity of 36 ± 0.9

Figure 10 Table 3 .Figure 9 .Figure 10 .

 103910 Figure 10 shows the distance 15 meters (a) and 12 meters (b) respectively. The attenuation efficiency of the 3 mm CdTe sensor in the energy region of 60 Co enables a timely localization of the test object with a maximum of 200 s in the case of 25 meters distance. It was possible to visually localize the isotope in the right arm of the structure at 20 meters.

of 1 .

 1 25 compared to a 2 mm CdTe sensor, 1.84 compared to a 1 mm CdTe sensor and up to 2.71 compared to a 0.45 mm CdTe sensor (Fig. 4, 137 Cs). In a scenario in which the possible source of radiation is unknown to the user, the extension of the field of view (FOV) of the camera by means of the possible distance between source and detector enables the first assessment of the situation in a timely matter. This depends on the attenuation efficiency and angular resolution of the detector. The integrated 3 mm CdTe detector measured 21.5 counts per second emitting from a 60 Co source with an activity of 2.8 ± 0.07 Gbq in 20-meters distance in an open environment. The acquisition time with visual identification by the user was 116 seconds.

Table 2 .

 2 Detector settings and configurations used for the measurements of 3.1 and 3.2.

	Sensor thickness	Bias (V)	Leakage current (μA)	Pixel pitch (μm)
	3 mm CdTe	-1200	2.5 ± 0.2	55
	2 mm CdTe	-600	8 ±1. 2	110
	1 mm CdTe	-450	2.1± 0.2	55
	0.45 mm CdTe	-230	43±4	55

Acknowledgments:

The authors would like to express their thanks to all collaborations partners of project "DURCHBLICK" (BMBF, grant no. 13N14328) in particular Andreas Ciossek of Telerob GmbH (Germany) for providing the robot and Frank Woodcock of ZIPPERMAST GmbH (Germany) for manufacturing the mounting solution of the iPIX. We would like to thank the Fraunhofer FKIE (Germany), the Bundesheer (Austria) and the organizers of the ENRICH for the opportunity of on-site testing of the system.

Funding: This research was funded by the German Federal Ministry of Education and Research (BMBF), grant number 13N14328.

Author Contributions: Conceptualization, J.F. and S.P.; methodology, J.F. and S.P.; software, J.F. and S.P.; validation, M.K.S., J.S.U. and J.F.; formal analysis, J.F.; investigation, J.F. and M.K.S.; resources, M.F., F.C., A.F. and V.S.; data curation, , M.K.S., J.S.U. and J.F.; writing-original draft preparation, J.F., S.P. and V.S.; writing-review and editing, J.F.; visualization, J.F.; supervision, S.P. and M.F.; project administration, S.P. and M.F.; funding acquisition, M.F.. All authors have read and agreed to the published version of the manuscript.

Conflicts of Interest:

The authors declare no conflict of interest.


Bingqian Liu

Côme Bissuel

Franc ¸ois Courtot

Céline Gicquel

Dominique Quadri

A hierarchical decomposition approach for the optimal design of a district cooling system

Keywords: Local energy system design, mixed-integer linear programming, Branch-and-bound, Hierarchical decomposition

A district cooling system is a centralized cooling supply system providing air conditioning to a set of buildings located in the same district. Designing and sizing such a system is very complex, as both the initial construction cost and the operation cost of the cooling system during its entire life must be considered. We first propose a modeling approach aiming at formulating this combinatorial optimization problem as a mixed-integer linear program of tractable size. We then extend a previously published hierarchical decomposition technique in order to find the optimal solution in an efficient way. Finally, we provide preliminary computational results based on a real-life case study located in China.

INTRODUCTION

A district cooling system (DCS) is a centralized cooling supply system. It consumes electricity to cool down water and distributes it through an underground pipe network to the buildings in the district to provide them with air conditioning. DCSs usually are highly energy-efficient cooling systems. Thus, according to the Electrical and Mechanical Services Department of Hong Kong (EMSD, 2020), using DCSs instead of traditional air-cooled air-conditioning systems results in energy savings of 35%. DCSs also compare well with individual water-cooled air-conditioning systems using cooling towers as the energy savings may be up to 1 20%. Furthermore, this lower energy consumption leads to lower greenhouse gas emissions, which helps reduce the environmental impact of the system.

Designing a DCS involves choosing the type and number of chillers to be installed as well as the ice storage capacity. These decisions should take into account not only the construction costs, but also the operation costs of the system during its whole lifetime. Computing these operation costs is a challenging problem. Namely, the demand for cooling power is highly variable and features a daily, weekly and yearly seasonality together with random variations. Moreover, due to technical reasons owing to the chillers, these operations costs are not at all proportional to the produced cooling power. Thus, in order to accurately estimate them, a detailed schedule describing, on a hourly basis, the on/off status and the load allocation of each chiller should be built for an horizon spanning a whole year. Furthermore, the deployment of a district cooling system is usually not a one-shot decision but rather a process in which investment decisions are made step by step, following the development of the district and the upward trend of the average demand over the years. This implies that a multi-phase strategic deployment plan should be built.

This optimization problem can be formulated as a mixed-integer program. However, its resolution poses several difficulties. The first one comes from the non-linearity of the chillers' performance curves. These performance curves give, for each chiller, the amount of electricity consumed as a function of the amount of produced cooling power and thus play a key role in the estimation of the system operation costs. Second, the need to simultaneously build a multi-year phasing plan and a detailed operational schedule for each day of the planning horizon leads to the formulation of a huge mathematical program, which cannot be solved directly by current mathematical programming solvers. Finally, the use of classical decomposition methods, such as the Benders' decomposition approach, is not straightforward as it would imply sub-problems involving binary and/or integer decision variables. The mixed-integer program modeling the problem is thus computationally intractable as such.

However, solution approaches exploiting the natural hierarchy between decisions relative to the system design and decisions relative to the daily operation schedules have been investigated in several works. [START_REF] Weber | Design and optimization of district energy systems[END_REF] thus studied the optimal design of a multi-energy system. A two-level optimization method is implemented. The master optimization level explores the set of possible system designs with the use of an evolutionary algorithm. For each considered system design, the slave optimization level calculates the optimal cost by linear programming. This approach does not provide a guaranteed optimal solution and the formulation of the slave sub-problems as linear programs does not allow to accurately model the way the system operates in practice. [START_REF] Iyer | Synthesis and operational planning of utility systems for multiperiod operation[END_REF] also use a bi-level method to optimize the choice and sizes of equipment for utility systems. At the design optimization level, the problem aims at fixing the system infrastructure. All binary decision variables relative to the operation schedules are removed from this master problem. Each time a potential infrastructure is found, the operation optimization level solves a set of single-period scheduling sub-problems taking the current system infrastructure as input data. Design cuts are used to tighten the gap between the solutions found at both levels. [START_REF] Yokoyama | Optimization of energy supply systems by MILP branch and bound method in consideration of hierarchical relationship between design and operation[END_REF] consider local energy supply systems and propose a customized Branch & Cut algorithm exploiting the hierarchical relationship between the design and operation decision variables of the mathematical program. The upper level problem corresponds to the initial optimization problem in which all operational integer and binary variables are kept but relaxed to be continuous. Each time an integer feasible solution is found at the upper level, a sequence of single-period independent operation sub-problems is solved to check the feasibility and value of the current design solution. Note that [START_REF] Iyer | Synthesis and operational planning of utility systems for multiperiod operation[END_REF] and [START_REF] Yokoyama | Optimization of energy supply systems by MILP branch and bound method in consideration of hierarchical relationship between design and operation[END_REF] both consider a single-phase variant of the problem, i.e. a variant in which all design decisions are one-shot decisions. Moreover, their design infrastructure did not allow short-term intra-day energy storage. This allows them to consider single-period (i.e. one-hour) scheduling sub-problems rather than multi-period (i.e. 24-hour) scheduling sub-problems, which significantly decreases the size of these sub-problems.

In the present work, we propose a solution approach for the optimal design, over a multi-phase investment horizon, of a local district cooling system in which intra-day ice storage is allowed. This approach relies on three key elements. First, we seek to reduce the size of the initial optimization problem. We thus consider a deployment plan involving a limited number of phases, some of which spanning several years. Moreover, we use the clustering approach described in [START_REF] Zatti | k-milp: A novel clustering approach to select typical and extreme days for multi-energy systems design optimization[END_REF] to select a small set of typical days to represent with the smallest possible loss of accuracy the various conditions under which the system will be operated. Second, we build a piecewise linear approximation of the performance curves of each chiller and propose a way to exploit their convexity to reduce the size of the formulation of the operation scheduling sub-problems. This results in the formulation of a large-size mixed-integer linear program (MILP). Thirdly, we develop an exact solution algorithm based on the hierarchical decomposition method recently proposed by [START_REF] Yokoyama | Optimization of energy supply systems by MILP branch and bound method in consideration of hierarchical relationship between design and operation[END_REF] to solve this MILP.

The remaining of the paper is organized as follows. Section 2 gives a detailed description of the problem under study. In Section 3, the approximation of the nonlinear performance curves and the clustering method used to identify typical days are first presented. The formulation of the optimization problem as an MILP is then provided. Section 4 introduces the hierarchical decomposition algorithm used to solve this MILP. In Section 5, the results of our preliminary computational experiments based on a real-life case study in China are reported.

PROBLEM PRESENTATION

Before the actual construction of a cooling system, a thorough and reliable system design is necessary. The system design consists in selecting the appropriate chillers and in sizing the ice storage capacity to be installed. The main objective is to design a system which will be able to meet the clients' cooling demand at all time while leading to the lowest long-term investments and operation costs.

Resources

A chiller is a machine that removes heat from a liquid by using a variety of techniques such as vapor-compression. We consider here electric-powered chillers which are used to cool water. These chillers can be classified into two main categories depending on their functionality. Standard chillers (denoted by STDC) only produce cooling power to satisfy the instantaneous demand of the customers. Ice chillers (denoted by ICEC) have two distinct operating modes: they either produce cooling power, but usually with a lower efficiency than the one of a STDC, or they produce ice. This ice can be stored for a few hours in an ice storage tank and be melted afterwards to provide cooling power. Within each category of chillers (STDC or ICEC), there are chillers with different production capacity levels. Each level corresponds to a predefined production range, i.e. to a minimum and maximum cooling power (or ice) it can provide per hour when turned on, and to a set of performance curves. Namely, the energy consumption of a chiller is a function of the produced cooling power and the ambient temperature. We thus have performance curves representing the relation between the electric power consumed by a chiller and the cooling power (or ice) produced under different ambient temperatures. Note that these performance curves are usually not linear. In the present paper, we will focus on the case of non-linear convex performance curves. The more general non-convex case is left for future work.

Another key resource in the system is the ice storage tank. This tank is linked to all the ice chillers of the system, can store the produced ice for a few hours and release it afterwards to produce cooling power. It has a maximal storage capacity, which can be chosen within a predefined range.

See Figure 1 for a schematic representation of the studied DCS.

Cooling demand

As highlighted in the introduction, the cooling demand to be satisfied by the local cooling system is highly variable. First, the demand for cooling power varies throughout the day and is usually much higher at daytime than at night. There are also weekly and yearly variations: the demand pattern of a weekday thus significantly differs from the one observed during the week-end and the total daily demand varies during the year, in particular with the summer and winter seasons. Finally, during the first operating years of the system, the demand displays a general upward trend as new clients join the district cooling system. After this initial period, the demand usually gets more or less stable for the rest of the system lifetime.

Parts of these demand variations display a seasonal pattern and are thus to some extent predictable. However, the demand is also subject to the presence of extreme weather conditions which are difficult to anticipate. Days with such effects should be considered separately from others. This translates in particular into the existence of extreme days, i.e. days in which the total daily or hourly demand is exceptionally low or high as compared to its usual value.

A key requirement is that the system should be able to satisfy the demand for cooling power at all time, whatever the hour of the day, the day in week or the season. In particular, it should be able to satisfy all the demand, even if it is exceptionally low or high.

Electricity supply

The chillers are powered by electricity bought from an external utility provider.

Similar to the cooling power demand, the electricity price displays daily, weekly and yearly variations. In particular, within a day, three types of periods, termed peak, flat and valley periods, can be distinguished. They correspond respectively to the highest, intermediate or lowest prices. Peak periods are usually at noon and in the evening, the valley periods at night and early in the morning and the rest of the day corresponds to flat periods. These price variations can be exploited to reduce the total energy cost, e.g. by producing ice at night when the demand for cooling power is low and the electricity rather cheap, storing it for a few hours and releasing ice to produce cooling power at daytime when the demand is high and the electricity more expensive.

Moreover, in many cases, the contract with the electricity provider includes a maximum allowed instantaneous power consumption from the grid. This upper limit is set contractually at the beginning of each year. The corresponding cost, called the contract fee, is proportional to the subscribed maximum power. This contract fee differs from the electricity consumption cost. It namely buys the permission of consuming electric power and sets an upper limit (expressed in kW) to this instantaneous consumption. In contrast, the electricity consumption cost depends on the total amount of consumed electric energy which is expressed in kWh.

Costs

The total cost of the system comprises two main parts: the design cost and the operation cost.

The design cost is the sum of the purchase, installation and maintenance costs of the chillers and ice storage tank, and of the annual contract fee. Once the cooling system design is chosen, this cost will be determined.

The operation cost corresponds to the cost of the electricity consumed when the system operates to satisfy the customers' cooling demand. Due to the nonlinearity of the chillers' performance curves, this cost is not proportional to the amount of produced cooling power. In order to accurately estimate it, we need to build a detailed production schedule for each day of the planning horizon using a fine time discretization.

The objective of this optimization problem is to find a system design which can minimize the sum of the design and operation costs. Therefore, the best system design as well as the most cost-efficient daily optimal operation strategy corresponding to the chosen system should be searched.

PROBLEM MODELING

Selection of typical days

The huge size of the optimization problem, which combines both long-term design and detailed (usually with an hourly timestep) production scheduling decisions, makes it computationally intractable. A possible way of getting over this difficulty consists in selecting, within the available data, a subset of typical days and extreme days which will represent the various conditions under which the system will be operated. These days should be carefully chosen as they will have a strong impact on the selection of the chillers and sizing of the ice storage tank.

This can be done by solving a clustering problem such as the k-medoid problem. In our case, this problem consists in partitioning the set of days belonging to the available historical time series into groups or clusters and in choosing one member (i.e. one day) in each cluster to represent it so as to minimize the total Euclidian distance between each day in the time series and its representative. In this work, we use the extension of the k-medoid problem proposed by [START_REF] Zatti | k-milp: A novel clustering approach to select typical and extreme days for multi-energy systems design optimization[END_REF] to select typical days for optimizing the design of our system.

Piecewise linear approximation of the performance curves

The performance curves of a given type of chiller are supplied by its manufacturer. They give, for a discrete set of values of the ambient temperature, the amount of electric energy consumed as a function of the output cooling power. Their non-linearity is a major source of difficulty for the resolution of our optimization problem. We thus propose to build a piecewise linear approximation of each curve.

This approximation is built using the following procedure. We first fix a predefined number of breakpoints in the piecewise linear approximation to be built. We then determine the coordinates of these breakpoints by heuristically solving a small non-linear optimization problem aiming at minimizing the distance between the approximate performance curve and the actual one.

Then, for each hour of each selected day, we use historical data about the average ambient temperature to determine the approximate performance curve that should be used to compute the energy consumption of each type of chiller during each time period.

Notation

The problem modeling relies on a three-level time discretization. The forecast lifetime of the system is first divided into a set of phases or investment periods, each one typically spanning one or several years: we assume that design decisions such as the installation of a new chiller can be only made at the beginning of a new phase. Let Φ be the number of considered phases. Within each considered phase, the various conditions under which the system will be operated will be represented by a number of preselected typical days and extreme days. We denote by D φ the set of typical days and extreme days used to represent the various daily demand patterns during phase φ ∈ {1, ..., Φ}. Each selected day d of phase φ has a weight w φ,d corresponding to the number of days it represents, i.e. to the number of days of the original historical time series which were assigned to the cluster it belongs to. Finally, in order to describe the intra-day variations of the cooling demand and electricity price, each selected day is divided into 24 one-hour periods. For the sake of readability, in what follows, we use the letter t to represent the time period (φ, d, h) corresponding to phase φ, selected day d and hour h. Let Dem t (resp. EP t ) represent the cooling power demand (resp. the electricity price) at time period t.

There are different types of chillers that may be installed in the system. Each type of chiller m = (p, l) can be described by its category p ∈ {ST DC, ICEC} (i.e. With respect to the design of the system, some restrictions have to be taken into account. There is namely an upper bound SD max m to the total number of chillers of type m that may be installed in the system. The total ice storage capacity built in the system should also stay below a maximum allowed capacity StoC max . Moreover, in terms of design costs, we denote by FC φ,m the fixed cost of investing in a chiller of type m at phase φ. This fixed cost comprises the construction or capital expenditure cost at phase φ and the total maintenance cost of the chiller over its whole lifetime (from phase φ to phase Φ ). Regarding the ice storage capacity, the installation cost is broadly proportional to the installed capacity. We denote by LC φ the cost of building one unit of ice storage capacity at phase φ. Finally, the unit subscription cost for the maximum instantaneous power is assumed denoted by SC φ . Note that all the design costs FC φ,m , LC φ and SC φ are assumed uniform along each deployment phase.

Variables

In order to model the problem as a mixed-integer linear program, we introduce two sets of decision variables.

The first set of decision variables corresponds to long-term design decisions which will determine the general structure of the system, together with the multiyear phasing plan. Thus, SD φ,m represents the integer number of chillers of type m to be installed at the beginning of phase φ, StoC φ is a continuous variable representing the ice storage capacity built at the beginning of phase φ and C φ is the maximum allowed electric power consumption contracted with the electricity provider for phase φ. These variables will be referred to as design decision variables in what follows.

The second set of decision variables are used to build the operational schedule for each selected day d ∈ D φ of each phase φ. In the present work, we will focus on the special case in which all performance curves of the chillers are convex.

This assumption allows us to use aggregate performance curves (see Appendix for more detail about this) and to build the schedule while considering, in each period t, aggregate scheduling variables, i.e. variables pertaining to the set of chillers of identical type producing the same commodity c during t, instead of disaggregate scheduling variables, i.e. variables pertaining to the status and output of each individual chiller installed in the system in period t.

We thus introduce S t,m,c the integer number of chillers of type m producing commodity c during period t, P t,m,c the total amount of commodity c produced by the chillers of type m during period t and Q t,m,c the total electric consumption of the chillers of type m producing commodity c in period t. Moreover, in order to monitor the ice storage and release, we define STO t the amount of ice stored in the tank at the beginning of period t and R t the amount of ice released during t. These decisions will be referred to as operation decision variables in what follows.

Objective function

We seek to minimize the sum of the design and operation costs of the system over its whole lifetime. The objective function of the mathematical program is thus given by: min

Φ ∑ φ=1 α φ ∑ m∈M FC φ,m SD φ,m + LC φ StoC φ + SC φ C φ + ∑ d∈D φ w φ,d 23 ∑ h=0 ∑ m∈M ∑ c∈C p EP φ,d,h Q φ,d,h,m,c (1)
where α φ is the actualization rate for the costs incurred in phase φ.

Constraints

Similar to the decision variables, the constraints of the mathematical model can be classified into two groups: design constraints and operation constraints.

Design constraints

Design constraints are constraints involving only design decision variables. In the present case, for each type of chiller m ∈ M , we have the following constraint which states that the total number of chillers of type m included in the DCS should be less than the maximum number of chillers of this type allowed.

Φ ∑ ϕ=1 SD ϕ,m ≤ SD max m (2)
Similarly, Constraint (3) makes sure that the total ice storage capacity of the system stays below the maximum allowed value. (3)

Operation constraints

Operations constraints are constraints involving operation decision variables, together with design variables in some cases. For each time period t, we have the following set of constraints.

Demand satisfaction

The demand for cooling power must be satisfied at all time. This can be done either by directly using the cooling power produced by the currently turned on chillers and by releasing some ice from the ice storage tank.

∑ m∈M P t, m,COLD + R t = Dem t (4)
Chillers Regarding the chillers, two sets of constraints should be introduced in the formulation.

First, for each type of chiller m ∈ M , the total number of operating (i.e.

turned on) chillers should be less than the number of chillers currently installed in the DCS:

∑ c∈C p S t,m,c ≤ φ ∑ ϕ=1 SD ϕ,m (5)
Second, for each type of chiller m ∈ M and each commodity c ∈ C p it can produce, the total amount produced by the turned on chillers should stay within the allowed production range: Moreover, the total amount of electric energy consumed in period t is limited by an upper bound, C φ × 1hour, which represents the maximum instantaneous power C φ contracted with the electricity provider times the duration of the period (one hour):

P t,m,c ≤ P max m,c S t,m,c (6)
P t,m,c ≥ P min m,c S t,m,c (7)
∑ m∈M ∑ c∈C p Q t,m,c ≤ C φ (9)
Ice storage Regarding the ice storage, three constraints should be considered in each time period. First, the amount of ice stored in the tank should not exceed the current storage capacity.

STO t ≤ φ ∑ ϕ=1 StoC ϕ (10)
Second, the amount of ice released during the period should be less than the amount stored in the tank at the beginning of the period.

R t ≤ STO t (11)
Third, the evolution of the ice inventory stored in the tank should comply with inventory balance equations. We first consider time periods (φ, d, h) corresponding to h ∈ {0...22} as the last hour of the day h = 23 requires a special treatment. For each period (φ, d, h) such that h ∈ {0...22}, we have:

STO φ,d,h + L ICEC ∑ l=1 P φ,d,h,ICEC,l,ICE -R φ,d,h = STO φ,d,h+1 (12)
Constraints (12) state that the amount of ice stored at the beginning of hour h + 1, STO φ,d,h+1 , is equal to the amount of ice already stored at the beginning of hour h, STO φ,d,h , plus the total amount of ice produced by the ice chillers of various capacity levels during h minus the ice melted during h. Note that the loss of ice stored in the tank during a day is assumed to be negligible.

Moreover, for each time period (φ, d, h) corresponding to h = 23, i.e. to the last hour of the selected day, we have the following inventory balance equation:

STO φ,d,23 + L ICEC ∑ l=1 P φ,d,23,ICEC,l,ICE -R φ,d,23 = STO φ,d,0 (13)
Namely, in practice, the entering inventory of a given day in the scheduling horizon is imposed by the leaving inventory of the previous day. In our case, we do not consider each individual day of the scheduling horizon but rather a number of representative days which will not necessarily occur successively in practice. We thus impose that the leaving ice inventory of a selected day d, computed as STO φ,d,23 + ∑ L ICEC l=1 P φ,d,23,ICEC,l,ICE -R φ,d,23 should be equal to the entering inventory of the same selected day d, STO φ,d,0 . This might be understood as the fact that the selected day d will be cyclically repeated w φ,d times in the simplified scheduling horizon used in our optimization problem for phase φ.

SOLUTION APPROACH

The mathematical program (1)-( 13) formulated in Section 3 displays a particular structure. We namely have a set of independent sub-problems. Each of them corresponds to optimizing the detailed schedule for a single selected day d of a single phase φ and involves operation variables and operation constraints relative only to the corresponding day (φ, d). These independent sub-problems are linked together by the design variables.

A Benders decomposition approach might seem appropriate for a problem displaying such a structure in which coupling variables link together a set of independent sub-problems. However, its application is not straightforward here as each sub-problem involves integer variables (namely variables S t,m,c ) and Benders decomposition algorithms rely on the strong duality theory to generate Benders cuts. We thus investigate another hierarchical decomposition approach recently proposed by [START_REF] Yokoyama | Optimization of energy supply systems by MILP branch and bound method in consideration of hierarchical relationship between design and operation[END_REF].

Hierarchical structure of the problem

To better highlight the hierarchical structure of the problem and more easily explain the decomposition approach, we will use a compact formulation of the problem.

Vector x represents in a synthetic way the design variables relative to all phases, i.e. x = (SD 1 , StoC 1 , C 1 , ..., SD φ , StoC φ , C φ ). Vector y k represents all the continuous operation variables relative to a given selected day k = (φ, d) whereas z k stands for all the binary or integer variables relative to day k.

With this notation, Problem ( 1)-( 13) can be formulated as follows.

min

Z = f 0 (x) + ∑ k∈K f k (y k , z k ) (14) h(x) ≤ 0 (15) g k (x, y k , z k ) ≤ 0 ∀k ∈ K (16) x ∈ Z ν (17) y k ∈ R µ ∀k ∈ K (18) z k ∈ Z λ ∀k ∈ K (19)
In the objective function ( 14), the term f 0 (x) corresponds to the design cost whereas the term f k (y k , z k ) computes the operation cost for each selected day k in the set K = {(φ, d), φ = 1...Φ, d ∈ D φ }. Constraints (15) correspond to the design constraints. Constraints ( 16) represent in a concise manner all the operations constraints relative to day k ∈ K . Constraints ( 17)-( 19) give the definition domain of each variable vector. Problem ( 14)-( 19) will be referred to as the Complete Model (CM) in what follows.

Note how all design variables are defined as integer variables in ( 17). This restriction is added to the problem as it is a necessary condition for the use of the hierarchical decomposition algorithm proposed by [START_REF] Yokoyama | Optimization of energy supply systems by MILP branch and bound method in consideration of hierarchical relationship between design and operation[END_REF].

This hierarchical decomposition algorithm uses as a starting point a relaxation of (CM) in which the design variables x are kept integer or binary whereas the binary or integer operation variables z are relaxed. This gives the following semirelaxed or upper level problem denoted by (SRM).

min

Z SRM = f 0 (x) + ∑ k∈K f k (y k , zk ) (20) h(x) ≤ 0 (21) g k (x, y k , zk ) ≤ 0 ∀k ∈ K (22) x ∈ Z ν (23) y k ∈ R µ ∀k ∈ K (24) zk ∈ R λ ∀k ∈ K (25)
Problem (SRM) thus involves the same number of variables and constraints as the initial problem ( 14)-( 19). However, the number of binary and integer variables is drastically reduced, which should ease its resolution.

Furthermore, the hierarchical decomposition algorithm relies on the key observation that in Problem (CM), when the design of the system is determined, the problem can be decomposed into a set of small independent operation subproblems. Let Problem OM k (x ) be the Operation Model relative to the optimization of the schedule of day k ∈ K , given a fixed design of the system described by vector x . It is formulated as:

min Z k (x ) = f k (y k , z k ) (26) g k (x , y k , z k ) ≤ 0 (27) y k ∈ R µ (28) z k ∈ Z λ (29)

Decomposition algorithm

The decomposition algorithm proposed by [START_REF] Yokoyama | Optimization of energy supply systems by MILP branch and bound method in consideration of hierarchical relationship between design and operation[END_REF] exploits the hierarchical structure described above. Thus, at the upper level, a relaxed version of the initial problem, i.e. problem SRM, is solved by a Branch & Cut algorithm. Each time a potential incumbent solution (x , ỹ, z) is found during this tree search, the corresponding values of the design variables x are used as input data to solve a series of independent scheduling sub-problems, namely problems OM k (x ), k ∈ K . This gives an accurate estimation of the feasibility and value of the potential design solution x at the operation level. If x is found to be feasible and less expensive than the current incumbent solution, it is accepted as the new incumbent solution. Otherwise, it is rejected. When all the branches are searched in the upper level Branch & Bound search tree, the current incumbent solution gives the optimal solution of the original problem. This hierarchical decomposition approach, provided it converges within the allotted computation time, thus guarantees the optimality of the found solution.

More precisely, the algorithm comprises the following main steps. See also Figure 2 for a flow chart of the decomposition algorithm.

Step 1 A Branch & Bound search tree is carried out within the feasible space of problem SRM: see 1 on the flow chart. The branching is done on the design variables x until there is no open node left, in which case we end the calculation, or an integer feasible solution of SRM, denoted by (x , ỹ, z), is obtained, in which case we go to Step 2. The value Z(x ) = f 0 (x )+∑ k∈K f k (ỹ k , zk ) provides a lower bound of the actual cost of the design solution x .

Step 2 Before accepting x as a new incumbent solution for the upper level problem, we check its feasibility and compute its actual cost Z(x ).

-We first initialize the value of Z(x ) as Z(x ) ←-Z(x ).

Then, for each k ∈ K :

-We solve Sub-problem OM k (x ) using a standard Branch & Cut algorithm: see 2 on the flow chart.

-If OM k (x ) is unfeasible, x cannot be feasible for Problem (CM). We stop and go to Step 4.

-Otherwise, we record the optimal integer solution of OM k (x ), (y k , z k ), and its optimal cost f k (y k , z k ).

-We update the current estimation of the actual cost of the design solution x by computing Z

(x ) ←-Z(x ) + f k (y k , z k ) -f k (ỹ k , zk ).
-If Z(x ) is larger than the incumbent value, x cannot be an optimal solution of Problem (CM). We stop and go to Step 4: see 3 on the flow chart.

-If all days in K have been considered, we go to Step 3. Otherwise, we go on with the next day in K .

Step 3 We replace the incumbent solution by (x , y, z) and the incumbent value by Z(x ) and go to Step 4: see 4 on the flow chart.

Step 4 We reject the current node in Problem SRM to prevent the solver from taking the semi-relaxed objective value of the design solution x , Z(x ), as a valid upper bound to be used for the rest of the Branch & Bound search at the upper level. This allows us to guarantee that the cutoff value taken into account to close nodes in the branching process is the actual cost Z(x ) of the current incumbent design solution x , see [START_REF] Yokoyama | Optimization of energy supply systems by MILP branch and bound method in consideration of hierarchical relationship between design and operation[END_REF]. Then go back to Step 1: see 5 on the flow chart. Finally, we noted in our preliminary experiments, that while running the hierarchical algorithm described above, we often had to solve multiple times the same operations problem OM k (.). Namely, two different design solutions x 1 and x 2 may e.g. be similar for the first phases and differ only for the last phases. It may also happen that two different design solutions x 1 and x 2 differ with respect to the decisions relative to the first phases of the horizon but give the same system infrastructure for the last phases. In both cases, operations problems OM k (x 1 ) and OM k (x 2 ) will be equivalent for all the selected days k corresponding to the phases in which x 1 and x 2 provide the same system infrastructure.

In order to avoid this useless computational effort, we thus modify the implementation of the hierarchical decomposition algorithm. When the operation cost of a newly encountered system infrastructure needs to be evaluated for a given phase, the operation sub-problems relative to this phase are solved and the corresponding optimal operation cost ∑ k st. d∈D φ f k (y k , z k ) is recorded in memory. Then, over the course of the algorithm, each time the operation cost of a design solution corresponding to the same system infrastructure in the same phase needs to be evaluated, we simply use the recorded value without recomputing it from scratch.

PRELIMINARY NUMERICAL RESULTS

In order to assess the performance of the proposed solution approach, we consider a real-life case study corresponding to a DCS under construction in China.

The expected lifetime of the cooling system is 30 years. The total annual cooling demand is anticipated to increase in the first three years and stay stable afterwards. We thus consider Φ = 3 investment phases. Phases 1 and 2 correspond to the first two years whereas Phase 3 corresponds to the last 28 years. To show the increase in demand, we compare in Table 1 the total yearly demand, which is the sum of demand of one year, and the maximum hourly demand for each phase. New chillers and storage capacity should be installed at the beginning of the first three years to guarantee that the cooling supply will meet the increasing demand. The hourly cooling demand value is predicted by combining historical data on the cooling consumption in the area and forecasts on the future number of clients which will connect to the DCS.

There are L ST DC = 3 types of standard chillers and L ICEC = 2 types of ice chillers available. Their maximum output capacity P max p,l,c and installation cost FC φ,m (in millions of CNY) are shown in Table 2. Note that chillers of type (ST DC, 1) and (ST DC, 2) have the same maximum output capacity P max p,l,COLD but chillers of type (ST DC, 2) are less efficient and less expensive than the ones of type (ST DC, 1). For each type of chiller and each corresponding commodity, the minimum output power P min p,l,c is equal to 0.10P max p,l,c . Figures 3 and4 display the performance curves of the available standard and ice chillers at an ambient temperature of 30 • C. We use B t,m,c = 4 breakpoints to build the piecewise linear approximation of the performance curve of the chillers of type m producing commodity c at time period t. The coordinates are determined by solving a small nonlinear optimization problem thanks to a heuristic method belonging to the numpy package in Python. Finally, for each type of chiller m, the maximum number of chillers that can be installed, SD max m is set to 10. The unit cost of installing ice storage capacity LC φ is 222.16CNY per kWh and the unit subscription cost for the maximum instantaneous power allowed, SC φ , is 276CNY per kW. Since the installed storage capacity and the subscribed contract power are discretized, a unit of storage capacity is 6000kWh and a unit of contract power is 3000kW. The installed storage capacity should be lower than StoC max = 200GWh and the contract power should be no more than 30GW. The discount rate α φ is 8%.

The electricity price features a daily seasonality but no weekly nor yearly variations. Figure 5 shows the unit price of electricity as a function of the hour of the day. As the electricity price does not change with the day in the week or the season, the typical days and extreme days are selected so as to represent as best as possible the variations in the cooling power demand. For each phase, 30 typical days are selected using the approach proposed by [START_REF] Zatti | k-milp: A novel clustering approach to select typical and extreme days for multi-energy systems design optimization[END_REF] and 4 extreme days are identified: the day with the highest hourly demand, the one with the highest total demand, the one with the lowest non-zero hourly demand and the one with the lowest non-zero total demand. Therefore, there are a total of 34 × 3 = 102 operation sub-problems.

This results in the formulation of a large-size MILP. The model involves 64381 variables, among which 22053 are integer and 600 are binary, as well as 101422 constraints. This MILP is solved using the hierarchical decomposition algorithm presented in Section 4. This algorithm is implemented in Python using the mathematical programming solver CPLEX12.8. All the problems are solved using a machine with a Intel Xeon 2.90GHz processor and 16GB RAM.

Table 3 describes the optimal solution in terms of the number of chillers and the storage capacity built in each phase, together with the contract power. Figure 6 shows the operation of chillers and storage in one of the typical days. The periods in which the ice production of (ICEC,2) seems to take negative value correspond to periods in which ice is produced and stored in the tank.

Furthermore, we carried out additional experiments in order to get a preliminary assessment of the computational efficiency of the hierarchical decomposition algorithm and to evaluate the impact of the selection of representative days on the infrastructure design. We thus created 4 additional instances based on our case study by varying the number of selected days per phase from 6 to 34. The corresponding results are provided in Table 4. The first line in the table indicates the number of selected days for each phase. Lines 2 to 5 provide indications on the size of the MILP to be solved: we note that the number of variables and constraints in the problem increases linearly with the number of selected days. Lines 6 to 9 correspond to the results obtained while directly solving Problem (CM) with CPLEX12.8 whereas Lines 10 to 13 correspond to the results obtained with the hierarchical decomposition algorithm. For each algorithm, we report the computation time (which was limited to a maximum of 2 hours), the total and design cost of the best solution found within the available time limit and the remaining gap (i.e. the relative difference between the values of the best solution and the best lower found after 2 hours of calculation).

These results first show that the hierarchical decomposition algorithm is significantly more efficient than the standard Branch & Cut algorithm embedded in CPLEX12.8. Namely, with this algorithm, we were able to solve to optimality the 5 considered instances, and this with a computation time divided by a factor of at least 2.25.

Moreover, results from Table 4 also show that the system infrastructure found by solving the MILP depends on the number of selected days. This can be seen among others by the fact that the investment cost varies with the number of selected days, which indicates that the structure of the system and the deployment plan vary with this number. We note however that, when using the hierarchical decomposition approach, both the total cost and the investment cost become stable when more than 26 selected days are considered for each phase. Additional computational experiments are needed to determine the exact value of the minimum number of selected days per phase above which the system design does not change any more.

Conclusion

We presented a modeling and solving approach for the optimal design of a district cooling system involving intra-day energy storage over a multi-phase horizon. The modeling approach relies on a clustering algorithm to identify a subset of typical days and on a piecewise linear approximation of the chillers' performance curves. It results in the formulation of a large-size mixed-integer linear program. An improved hierarchical decomposition method is then implemented to optimally solve this MILP. This decomposition method exploits the hierarchy between upper-level infrastructure design decision and lower-level operation scheduling decisions. Our preliminary computational results carried out on a real-life case study located in China show that the proposed hierarchical decomposition algorithm significantly outperforms a mathematical programming solver at providing optimal solutions of the MILP. Moreover, our results also show that, provided a minimum number of typical days are taken into account to estimate the operation costs, the final system infrastructure and the deployment phasing do not change with the subset of selected typical days, which is an important point to gain the trust of the decision makers.

There are several possible directions for future research suggested by the present work. First, additional computational experiments are needed to evaluate the impact of the use of a limited number of representative days and of the piecewise linear approximation of the chillers' performance curves on the solution provided by the proposed approach. Second, on a longer term perspective, it would be interesting to extend the proposed approach to consider non-convex performance curves for the chillers and to study more general local energy systems, in particular systems simultaneously providing a district with several sources of energy (electricity, heat, cold, ...).

Figure 1 :

 1 Figure 1: Overview of a DCS

  standard or ice chiller), which defines the list of commodities C p ⊂ {COLD, ICE} it can produce, and its production capacity level l ∈ {1, ..., L p }. Let M = {(p, l)|p ∈ {ST DC, ICEC}, l ∈ {1, ..., L p }} be the set of all chiller types. We denote by P min m,c and P max m,c the minimum and maximum output power of a chiller of type m = (p, l) producing commodity c ∈ C p . As explained above, the performance curve of a chiller of type m producing commodity c at time period t is approximately represented by a piecewise linear function comprising B t,m,c breakpoints. Let a t,m,c,b and o t,m,c,b be the abscissa and ordinate of breakpoint b of this piecewise linear function.

  Electric consumption As shown in Appendix, when the chillers' performance curve are convex, the aggregate electric consumption of the set of chillers corresponding to a given type m = (p, l) ∈ M and producing commodity c ∈ C p in t belongs to the epigraph of a piecewise linear function defined by breakpoints b = 1...B t,m,c . We thus have, for each m = (p, l) ∈ M , each c ∈ C p and each b ∈ {1...B t,m,c -1}, the following inequality: Q t,m,c ≥ s t,m,c,b P t, m, c + c t,m,c,b S t,m,c (8) where s t,m,c,b = o t,m,c,b+1 -o t,m,c,b a t,m,c,b+1 -a t,m,c,b is the slope of the line segment between breakpoints b and b + 1 and c t,m,c,b = o t,m,c,bs t,m,c,b a t,m,c,b the corresponding constant value. Note that Constraints (8) accurately compute the electric consumption of a set of chillers only if the corresponding performance curve is convex. The more general case of non-convex performance curve is left for future work.

Figure 2 :

 2 Figure 2: Flow chart of hierarchical decomposition algorithm

Figure 3 :

 3 Figure 3: Performance curves of STDC at 30 • C

Figure 5 :

 5 Figure 5: Electricity Price of a Day

Table 1 :

 1 Phases and Demand

	20

Table 3 :

 3 Optimal system design and phasing obtained with 30 typical days and 4 extreme days at the operation level

Table 4 :

 4 10 8 5.70 × 10 8 5.75 × 10 8 5.75 × 10 8 5.75 × 10 8 Invest. Cost 2.45 × 10 8 2.35 × 10 8 2.36 × 10 8 2.36 × 10 8 2.37 × 10 8 Value 5.50 × 10 8 5.70 × 10 8 5.74 × 10 8 5.75 × 10 8 5.75 × 10 8 Invest. Cost 2.45 × 10 8 2.35 × 10 8 2.35 × 10 8 2.35 × 10 8 2.35 × 10 8 Comparison between complete model and hierarchical decomposition

	Instance	#Days	6	10	18	26	34
	Size	#Var	11461	19021	34141	49216	64381
		#IntVar	3909	6501	11685	16869	22053
		#BinVar	180	240	360	480	600
		#Const.	17938	29854	53710	77566	101422
	CPLEX	Comp. time	1726s	>7200s	>7200s	>7200s	>7200s
		Obj. Value 5.50 × Relative Gap 0.00%	0.02%	0.15%	0.34%	0.42%
	Hier. Dec. Comp. time	220s	922s	942s	1435s	1831s
		Objec. Relative Gap	0.00%	0.00%	0.00%	0.00%	0.00%

APPENDIX Convex performance curves

When the performance curves of the chillers are convex, we have the following property.

For a set of identical chillers producing the same commodity, the optimal load allocation consists in equally distributing the total output power between the chillers.

Proof. The proof is done by contradiction. Suppose we have two identical chillers, producing a total cooling power of P with chiller 1 producing P 1 and chiller 2 producing P 2 > P 1 . Let π : P → Q = π(P) be the convex performance curve of these two chillers. The total amount of electricity consumed by the two chillers producing P is Q = π(P 1 ) + π(P 2 ).

We show that this load allocation is not optimal, i.e. that it is possible to reduce the total amount of consumed electricity. Namely, let δP be a small variation in the output. By decreasing the output of chiller 2 by δP, we can obtain a decrease in the electricity consumed by this chiller of π (P 2 )δP, where π is the derivative function of π. In order to still be able to provide a total output of P, we increase the output of chiller 1 by δP, which leads to an increase in its electricity consumption of π (P 1 )δP.

By convexity of function f , π (P 1 ) < π (P 2 ). Hence the total amount of electricity consumed with the load allocation (P 1 + δP, P 2 -δP) is smaller than the one consumed with the load allocation (P 1 , P 2 ). The result follows.

Let us now focus on the case in which the performance curve π is convex and piecewise linear. When multiple identical chillers are simultaneously producing the same commodity, the relation providing the total amount of consumed electricity as a function of the total amount of output power can be plotted as an aggregate performance curve. We have the following property:

Let π be the piecewise linear and convex performance curve of a given type of chiller. π involves B breakpoints. Let (a b , o b ) be the abscissa and ordinate of breakpoint b.

The aggregate performance curve Π S of S identical chillers of this type is also piecewise linear and convex. It involves B breakpoints whose coordinates are given by (Sa b , So b ).

Proof. Let us consider the case where the total output of the S chillers is P ∈ [SP min ; SP max ] where [P min , P max ] is the production range of a single chiller.

By Lemma 1, the optimal load allocation consists in requiring each chiller γ = 1...S to produce the same output P γ = P S . Let b be the index of the breakpoint of function f such that P γ ∈ [a b , a b+1 ]. The energy consumed by each chiller γ is thus given by: Q