Applied Optimal design of multi-energy systems with seasonal storage

(cid:129) Novel MILP approaches to enable design of MES including seasonal energy storage. (cid:129) Good accuracy and much lower computational complexity compared to current approaches. (cid:129) Realistic Swiss case-study evaluated in terms of total annual cost and emissions. (cid:129) Extensive sensitivity analysis de ﬁ ning design guidelines for seasonal energy storage. Optimal design and operation of multi-energy systems involving seasonal energy storage are often hindered by the complexity of the optimization problem. Indeed, the description of seasonal cycles requires a year-long time horizon, while the system operation calls for hourly resolution; this turns into a large number of decision variables, including binary variables, when large systems are analyzed. This work presents novel mixed integer linear program methodologies that allow considering a year time horizon with hour resolution while sig- ni ﬁ cantly reducing the complexity of the optimization problem. First, the validity of the proposed techniques is tested by considering a simple system that can be solved in a reasonable computational time without resorting to design days. Findings show that the results of the proposed approaches are in good agreement with the full-scale optimization, thus allowing to correctly size the energy storage and to operate the system with a long-term policy, while signi ﬁ cantly simplifying the optimization problem. Furthermore, the developed methodology is adopted to design a multi-energy system based on a neighborhood in Zurich, Switzerland, which is optimized in terms of total annual costs and carbon dioxide emissions. Finally the system behavior is revealed by performing a sensitivity analysis on di ﬀ erent features of the energy system and by looking at the topology of the energy hub along the Pareto sets.


Introduction
Recently, the energy sector has been riding a wave of grand transformation: the necessity of decreasing the environmental impact has led to the deployment of conversion and storage technologies based on renewable energy sources [1]. In this context, multi-energy systems (MES) represent a new paradigm which exploits the interaction between various energy carriers (e.g. electricity and heat) at design and operation phase, allowing for improved technical, economic and environmental performance of the system [2]. Within this framework, seasonal storage systems have recently caught much attention due to their ability to compensate the seasonal intermittency of renewable energy sources [3]. However, compensating renewables fluctuations at the seasonal scale is particularly challenging: on the one hand, a few systems, such as hydrogen storage and large thermal storage, allow offsetting seasonal variations in renewable energy generation; on the other hand, the optimal design and operation is complicated by the large number of decision variables, due to the required length and resolution of the time horizon.
Several works provide comprehensive reviews of the model formulations and computer tools adopted for investigating MES and their integration with renewable energy sources and storage technologies. For instance, Alarcon-Rodriguez et al. focused on the multi-objective planning of distributed energy resources [4]; Connolly et al. presented a review of the computer tools implemented for analyzing the integration of renewable energy into various energy systems [5], whereas Keirstead et al. [6] and Allegrini et al. [7] focused on urban energy system models; Mancarella provided an overview of concepts and models for the planning and analysis of multi-energy systems [2]. When storage technologies are available, the optimal design of MES is significantly complicated by the necessity to consider the system operation already at design phase to accurately make use of the storage systems. Although a few nonlinear approaches have been proposed, for instance by Elsido et al. [8], mixed-integer linear programming (MILP) has been particularly favored as optimization framework for MES design and operation since it catches well the features of the system with a reasonable computational complexity. The problem of optimal technology selection and unit commitment through MILP formulation has been extensively investigated in the past. For example, Marnay et al. presented the case of a commercial building micro-grid with heat and electrical storage [9]; Hawkes and Leach extended the study by considering a hospital and residential buildings [10]. Later, Angrisani et al. investigated the energy, economic, and environmental performance of micro tri-generation systems [11]. Fazlollahi et al. introduced methods for multi-objective design of complex energy systems [12], and Ahmadi et al. presented the thermodynamic modeling and multi-objective optimization of an energy system for the simultaneous generation of electricity, heating, cooling and hot water [13]. Whereas these works were mainly focused on small, yet centralized systems (i.e. one hub for different end users), a number of studies also investigated energy distribution among the different nodes of decentralized energy systems (i.e. multiple hubs for different end users). For example, Genon et al. Greek letters first principle-to-electrical efficiency ratio ( presented the environmental assessment of small district heating systems [14], while Weber and Shah optimized the structure of the heating network in addition to technology selection and unit dispatch [15]. Furthermore, various case-studies have been proposed by applying such tools: Mehleri et al. investigated the optimal design and operation of different neighborhoods within the Greek energy sector [16,17]; Omu et al. focused on the possible measures to reduce the carbon footprint of the UK energy sector [18], and Bracco et al. presented the design and operation of the micro-grid built at the University campus in Savona, Italy [19]. Finally, the interest in multi-energy systems at different levels (from residential to national) has triggered the development of commercial tools for MILP formulation. These include EnergyPlan, an optimization tool developed at Aalborg University, which simulates the operation of national energy systems [20], and DER-CAM, developed by the Lawrence Berkeley National Laboratory [21]. Recent improvements of DER-CAM tool include building retrofitting [22] and detailed models of sensible thermal storage systems [23]. In these works, the remarkable complexity of the optimization problem required significant model simplifications. In specific, all these past studies considered a one-year time horizon based on design days. Recently, Elsido et al. investigated the weekly periodicity in the operation of CHP units with thermal storage, describing the entire year through typical weeks [8]. These approaches are suitable when analyzing intraday storage installation and operation, whilst they are not able to describe seasonal operation cycles due to the discontinuity between the selected design days (or weeks). Longer time horizons have been considered when investigating only the operation of multi-energy systems, e.g. in [24], but not at design phase.
This work aims at addressing and overcoming this issue. To this end, two novel MILP formulations based on the coupling of typical design days and on the distinction between two groups of decision variables are proposed. These approaches allow considering an hour resolution of the entire year while maintaining a low number of binary variables, thus enabling the solution of complex systems featuring seasonal energy storage. In the following, first, the proposed techniques are tested by considering a simple system that can be solved in a reasonable computational time without resorting to design days (i.e. year time horizon with full hourly resolution). Next, the presented methodology is applied to a residential multi-energy system based on a neighborhood in Zurich, Switzerland: different hub designs are determined and analyzed to minimize the total annual cost and CO 2 emissions. Finally, the operation of the energy storage during the whole year is presented and discussed.
The paper is structured as follows. Section 2 describes the features of the investigated multi-energy system, formulates the traditional optimization problem and describes the novel approaches. Section 3 presents the results. Finally, in Section 4 conclusions are drawn.

System description and formulation of the optimization problem
The multi-energy system considered in this study has the primary objective of supplying the energy demanded by a specified user. The MES is connected to the gas and electrical grids and includes a set of conversion technologies, both traditional and renewable-based, and of storage units. Fig. 1 provides a schematic representation of such a system. The weather conditions, electricity and gas prices, and energy demand profiles are the three sets of inputs to the optimization problem. The model returns the optimal system design in terms of technology selection and size, and the optimal operation of the installed units.
The problem is mathematically formulated as a MILP, where binary variables are introduced to model the performance of the conversion units, and the capital cost of the investigated technologies. The MILP can be written in general form as where c and d represent the cost vectors associated to continuous and binary decision variables, x and y respectively; A and B are the corresponding constraint matrices and b is the constraint known-term; N x and N y indicate the dimension of x and y, respectively. In the following, the optimization problem is described in detail in terms of input data, decision variables, constraints, and objective function.

Input data
The input data are time-dependent profiles for 2016 in a fraction of Altstetten, a neighborhood within the city of Zurich, Switzerland, which features a peak electricity demand of 0.43 MW and a peak thermal demand of 2 MW. The hour profiles for all input data are reported in the appendix. Data are assumed to be constant along the lifetime of the system and known without uncertainty, i.e. it is assumed that the possible realizations of the uncertainty and the evolution of the demand profiles are represented by the historical data. Indeed, the optimization of MES considering the uncertainties of the input data (e.g. electricity price, heat and electricity demand, renewable potential) is an important and open topic, which is being tackled by several researchers by means of different methodologies (e.g. [25][26][27][28][29][30][31][32]). However, such research activity is not within the scope of this paper, which aims at providing a MILP framework for implementing seasonal storage, regardless of the extent to which input data are representative of the future. Moreover, findings of this work can be extended and implemented in optimization problems under uncertainties. In this framework, solving the optimization problem at every time instant of the time horizon, i.e. every hour of the year, provides the optimal solution to the problem. In our MILP implementation, inputs to the optimization problem are: i. The expected weather conditions, namely air temperature,  ∈ Θ T , and solar radiation, where T indicates the length of the time horizon. Data are available at every hour of the year, hence T = 8760. ii. The expected prices of energy utilities, namely import prices, u e and  ∈ u T g , for electricity and gas, respectively, and export electricity price, , variable for different technologies, where M indicates the number of available technologies. While the import gas price is considered to be fixed along the year, the import electricity price is not constant. Notably, the electricity price features a very volatile profile, which is the result of a complex market regulation together with a broad portfolio of electricity suppliers (i.e. renewables-, nuclear-, fossil fuel-based). The Swiss electricity spot market price for 2016 with an hourly resolution is considered. iii. The expected electricity,  ∈ L T e , and heat,  ∈ L T h , demand profiles. Data are available at every hour of the year. iv. The set of available technologies with the corresponding performance and cost coefficients.

Decision variables
The following decision variables are returned: i. The size of the installed technologies,  ∈ S M (deciding S also implies selecting the technologies). ii. The on/off status, ∈ . Note that the exported electrical power must be determined separately for each technology due to the different selling prices.
While the size and selection of the installed technologies is typically referred to as a design variable, decision variables (ii)-(v) are denoted as operation variables and are determined at every hour of the year.

Constraints
The constraints of the optimization problem can be divided into two categories: (i) performance of the conversion and storage technologies, and (ii) MES energy balances. Both are described in detail in the following.
Performance of conversion and storage technologies. The set of available technologies is indicated with M and includes all the technologies described below. Here, affine and piecewise affine (PWA) correlations describing the performance of the considered technologies are derived based either on first-principle models or on manufacturer data. In particular, the approach proposed by Yokoyama et al. [33] is used to linearize the power and size dependency of conversion performance. In the following equations, S refers to the size of the technology of interest, i.e. the rated input power. For each technology M ∈ i , S i must vary between a minimum and a maximum value, S i min and S i max respectively: where the binary variable ∈ a {0,1} i indicates whether the i-th technologies is installed. The following constraints hold for all time steps ∈ … t T {1, , } (the linear coefficients for all the investigated technologies are reported in Table 1): • Photovoltaic (PV) panels, generating electricity from solar energy: the generated electrical power P t is expressed as where S is the installed PV area, limited by the roof area availability; η is the pre-determined conversion efficiency, modeled as a function of air temperature Θ t and solar radiation I t , based on the approach suggested by De Soto et al. [34]: an average value around 0.15 is obtained for our case using that approach. Note that P t is the t-th element of the vector P (similar notation is used for all the other vectors).
• Thermal solar (TS) panels, generating heat from solar energy: the generated thermal power P t is expressed as where S is the installed TS area, limited by the roof area availability; η is the conversion efficiency, which is typically given as function of the generated heat temperature: here it is assumed to be constant, and equal to 0.65 [35]. Note that the sum of installed PV and TS area is also limited by the total availability.
• Electricity-driven heat pumps (edHPs), generating heat by absorbing electricity: given the Swiss boundary conditions, here the heat pumps are only used to provide heat. The generated thermal power is expressed as where Here x t is a binary variable indicating whether the device is turned on at time step t, producing power but also incurring in a minimum power consumption + δS ζ ; the parameters α β γ δ ζ κ ν , , , , , , linearly correlate the generated power P t to the input power F t and the unit size S, and are determined by fitting manufacturer data for several sizes and part-load conditions [36]; the influence of air temperature on the efficiency is taken into account through the given function f (Θ ) t , based on the data presented in [37]. Finally, note that the bilinearity Sx t between a continuous and a binary variable can be formulated as a linear constraint through the introduction of an auxiliary variable = ∼ S Sx t t [38]: • Micro-gas turbines (MGTs) generating electricity and heat from natural gas (NG): the generated electrical power is expressed as where Here two different sets of linear coefficients αν are implemented to distinguish between devices below and above a certain size S int . In other words, different affine approximations are implemented for different size ranges, which are identified through the binary variable ∈ b {0,1} 2 . This allows a better fitting of manufacturer data [43]. Micro gas turbines produce at the same time electrical and thermal power Q t , the two being related by where ρ is the average ratio of first-principle efficiency to electrical efficiency.

• Fuel cells and electrolyzers: both solid oxide fuel cells (SOFCs) and
proton exchange membrane fuel cells (PEMFCs) are considered. While SOFC run on NG, PEMFC can use NG or hydrogen (H 2 ) as fuel.
A proton exchange membrane (PEM) electrolyzer is also included, generating hydrogen and oxygen by absorbing electrical power. The produced hydrogen and oxygen can be stored in pressurized tanks. The possibility of producing and re-converting hydrogen is here referred to as power-to-gas (PtG). For the fuel cell, P t and F t refer to generated electrical power and inlet fuel power (fuel LHV), respectively. For the electrolyzer P t and F t refer to generated fuel power (hydrogen LHV) and absorbed electrical power, respectively.
For both devices, a piecewise affine approximation based on firstprinciple models is implemented to describe the conversion efficiency, as discussed in [44,45]. For all line segments = … i m {1, , }: Here the parameters α β and i i are the coefficients of the i-th approximating line segment. Note that due to the modularity of these systems the size does not significantly affect the performance, translating into = = = γ ζ ν 0. Like MGTs, fuel cells produce at the same time electrical and thermal power Q t , the two being related by Eq. (13). Note that startup/shut-down and ramp-up/ramp-down limitations of fuel cells, electrolyzers and micro gas turbines are modeled through the set of linear constraints presented by Arroyo and Conejo in [46] (in particular, we refer to Eqs. (1)- (12) in that paper), where two additional binary variables y t and z t are introduced. An additional remark about the oxidant fed to the fuel cells (FCs) is worth making. In this work, oxygen from air is used instead of pure oxygen produced by the electrolyzer. Although this translates into a lower electrical efficiency of the fuel cell, it avoids: (i) the oversizing of electrolyzer and storage tanks (oxygen is produced by the PEM electrolyzer with a H 2 :O 2 ratio of 2:1 and consumed by the PEM fuel cell with a H 2 :O 2 ratio of 1.15:1), (ii) the injection of excess hydrogen into the natural gas grid, which is still a debated procedure in terms of cost and flow rate limitations.
• Boiler, generating heat from natural gas: the generated thermal power is expressed as Here, the size dependency of boiler performance is neglected, and no minimum value is imposed for the generated power. A conservative treatment is adopted, as the boiler is used as a reference technology for thermal generation.
• Storage systems: three types of storage systems are considered: (i) hot water sensible thermal storage (HWTS), (ii) lithium battery (LiB), and (iii) H 2 tank storage (HS). Although latent and thermochemical thermal storage systems promise to be suitable solutions for long-term storage in the near future [47], sensible thermal storage represents the cheapest, most developed and commercial system. In this work, HWTS is considered, where the water is stored at 90°C and cooled to 65°C. Due to the fairly high energy losses and the low energy density, this system is mainly used to overcome short-term mismatch between thermal energy generation and use. Similarly, battery storage is mainly used for short-term electricity compensation due to its energy losses and high investment cost. On the contrary, H 2 storage features negligible energy losses, and therefore represents a promising solution for offsetting seasonal mismatch between renewable energy generation and energy consumption. The HS is coupled with the PtG system to generate, store and use hydrogen. In this context, hydrogen can be generated at high pressure, 40 bar here. While on the one hand HS has very limited energy losses in time, on the other hand the round-trip efficiency of the overall PtG process is around 40%, much lower than that of thermal and battery storage systems (see Table 1). For this reason, HS represents an efficient alternative when storing energy for long time periods. It is worth mentioning that the different utilization of the storage units typically translates into a larger installed size for the HS than for HWTS and LiB. All the aforementioned storage units are described through the following linear equations: Here P t assumes positive or negative values for charging and discharging power, respectively; Λ and Π are self-discharge parameters, and g(Θ t ) expresses the influence of the ambient temperature on the losses of the storage devices, as suggested in [23] for the thermal storage; η indicates the charging/discharging efficiency; t Δ is the duration of the time interval t τ ; is the time required to fully charge or discharge the storage. Note that the same value for charging and discharging efficiency is assumed. As the two efficiencies are typically similar, there is no need for binary variables specifying whether the device is charging or discharging. The periodicity constraint, Eq. (19), imposes the same storage level at the beginning and at the end of the year.
Note that in the present form the optimization framework does not include considerations on possible malfunction of the devices.
MES energy balances. The set of considered energy carriers within the multi-energy system is indicated with N and includes: • Electricity (e): it can be consumed and generated by the conversion technologies, stored by the storage units, exchanged with the electrical grid and utilized by the end users.
• Heat (h): it can be generated by the conversion technologies, stored by the storage units and utilized by the end users.
• Natural gas (g): it can be imported from the gas grid and consumed by the conversion technologies.
• Hydrogen (H 2 ): it can be consumed and generated by the conversion technologies, stored by the storage units, and injected into the gas grid.
It is worth noting that cooling is not considered here; however, it could easily be included.
The sum of imported and generated power must equal the sum of exported and used power for each energy carrier N ∈ j for all time steps ∈ … t T {1, , }. This can be expressed in general form as Here, U is the imported energy, P the generated energy, V the exported energy, F the absorbed energy and L the energy required by the end users. With this notation, the set of technologies M also includes the direct utilization of the energy carriers. Note that limitations on the imported and exported energy could easily be included if necessary. An additional remark about the modeling of electricity and heating demands is needed. These demands are modeled for several buildings based on a neighborhood in Zurich, Switzerland, following the approach presented in [48]. A centralized energy hub is assumed to be located at the center of the neighborhood, and an approximate grid topology for the district heating is considered using a shortest path algorithm. Based on the calculated distance between the hub and the buildings, the original heat demand of each building is updated to account for transport losses. A specific heat demand increment of 4% per km is assumed [49]. Also, an auxiliary electricity consumption of 8.5% of the heat demand is considered for pumping [15]. Finally, the electricity and heat demands for all buildings are aggregated and the total demands are implemented.

Objective function
The objective function of the optimization problem is the total annual cost of the system, J, given by the sum of three contributions, namely the capital, J c , operation, J o , and maintenance, J m , annual costs.
The annual capital cost is expressed as where λ i and μ i represent the variable and fixed cost coefficients for the i-th technology; S i indicates the unit size, i.e. the rated input power for non-solar conversion technologies, the rated stored energy for storage technologies, and the installed area for solar panels. To compute the equivalent annual investment cost, the annuity factor ω is included, which is calculated using the standard definition and based on an interest rate of 6%. A piecewise affine approximation is implemented to model the size dependency of the capital cost. For each technology M ∈ i : are the minimum and maximum size values of every line segment, respectively. Since the capital cost must be minimized, and due to the concavity of the curve, a binary variable for each line segment is required to identify the active linear approximation. The cost approximation for two exemplary technologies, namely heat pump and thermal storage, is illustrated in Fig. 2. Table 2 reports the cost coefficients for all the investigated conversion and storage technologies, as well as their capacity constraints.
The annual operation cost is calculated based on the amount of imported and exported electricity and gas during the year: where import and export prices, u and v, and powers, U and V, depend €/kWh ∀ t, while a lower electricity selling price, 0.05 €/kWh, is assumed for the other conversion and storage technologies. The Swiss electricity spot market price for 2016, which varies from 0.01 to 0.12 €/kWh is used [56], whilst a constant import price of 0.064 €/kWh is used for natural gas [57].
The annual maintenance cost is given as a fraction ψ of the annual capital cost [52,53]: Along with the total annual cost, the system is also evaluated in terms of environmental performance by means of the ε-constraint method. This translates into a minimum-cost optimization problem, where the annual CO 2 emissions are constrained below a maximum threshold. The total annual emission is calculated as the sum of the contributions of electricity and gas imported from the grid: where ε j is the specific carbon dioxide emission of carrier j, assumed to be constant along the year. The values of the emission coefficients for electricity and gas are = ε 0.137 e ton CO2 /MWh and = ε 0.237 g ton CO2 /MWh, respectively. Note that the Swiss energy mix strongly relies on hydroelectric and nuclear power, thus resulting into low CO 2 emissions for grid electricity.

Time horizon modeling
So far, we have described the optimization problem considering as time horizon one year with hour resolution, i.e. = T 8760. Because of the large number of variables and constraints arising when investigating different technology options, the full scale optimization is hardly feasible and the time variability during the year is normally simplified. Recently, Pfenninger presented a systematic analysis of the different techniques to reduce the time resolution of energy models [58], including time slices and representative design days or weeks. Here, a set of D design days and their sequence σ along the year are identified by using the MATLAB®-embedded clustering algorithm kmeans [59]. Based on weather conditions, energy prices and load profiles (i.e. input data (i)-(iii)), the algorithm solves a mixed-integer nonlinear problem to identify the most representative set of D design days and to assign every day of the year to a specific design day ∈ … d D {1, , }, characterized by its own typical hour resolution. The original formulation of k-means clustering, or Lloyd's algorithm, was presented by Lloyd in [60]. It is worth mentioning that the clustering procedure would level off the original input profiles, reducing the variation between the minimum and maximum values. Consequently, a system design based on such clustered input profiles would not be able to meet the energy demand at every hour of the year. To address this issue, and having in mind that the primary goal of the hub is to satisfy the user demand, the clustering procedure has been constrained to maintain the minimum and maximum values of the original demand profiles.
The traditional MILP formulations based on design days are not able to include seasonal storage. In the following, three methods for modeling the time horizon are discussed: M0, which is the existing method that excludes seasonal storage, and M1 and M2, which are proposed in this work as alternative solutions that allow including seasonal storage.
Method 0 (M0) -Traditional approach. The traditional approach, hereafter denoted as method 0 (M0), considers uncoupled design days: The sequence of design days along the year is not considered in the optimization problem, which is therefore formulated independently for each design day d. As an example, the equations describing the power generated by fuel cell and electrolyzer devices (Eqs. 14,15) are recast as: where ∈ … d D {1, , } indicates the d-th design day and ∈ … k K {1, , } the k-th daily time instant, with hour resolution hence = K 24. Similarly, the storage dynamics (Eqs. (18) and (19)) becomes where the periodicity constraint, Eq. (33), is written independently for each d-th design day, thus decoupling each day from both the previous and the next. Note that this may result into a mismatch between the energy stored at the end of a day and the energy stored at the beginning of the following one. Following the same logic, the energy balances of Eq. (23) are recast as: and the annual operation cost of Eq. (27) is given by ). Although this approach reduces the complexity of the optimization problem, as D is smaller than 365, the discontinuity between different design days does not allow to describe seasonal operating cycles, therefore resulting in a daily operation of the multi-energy system. In the following, two novel approaches are proposed to allow for longterm system operation while at the same time limiting the computational complexity of the optimization problem.
Method 1 (M1) -Coupling design days. The first approach proposed in this paper, hereafter named method 1 (M1), allows for the optimization of the level of stored energy hour-by-hour while still describing the year time variability through typical design days. This requires the introduction of the sequence σ of design days, which allows to couple successive days of the year. In other words, the energy stored during the last hour of every day of the year is connected to the energy stored during the first hour of the following day, where each day of the year is described by a given design day. Note that σ takes each day of the year y and returns the corresponding design day d, i.e. ⩽ ⩽ ⩽ ⩽ y σy D 1 365, 1 ( ) . In this case, the storage dynamics of Eqs. (32) and (33) is written as: where ∈ … y Y {1, , } indicates the y-th day of the year, with = Y P 365; σ y k ( ), is the charging/discharging power at day y, described by the typical day σ y ( ). Eqs. (36)- (38) imply that, although two days of the year described by the same design day must be characterized by the same variation in stored energy, they can have different values of stored energy at the beginning of each day. The improved storage flexibility obtained with respect to M0 comes at the price of a larger number of variables and constraints. Indeed, D K ) can now assume different values at every hour of the year, and the corresponding number of constraints is imposed by Eqs. (36) and (37). As for the full scale optimization, the periodicity constraint, Eq. (38), is imposed on the whole year, instead of on each design day.
Method 2 (M2) -Detailed input data. The second approach proposed, hereafter named method 2 (M2), allows to optimize the storage operation hour-by-hour while considering the actual energy demands.
The underlying consideration is that the computational complexity of the optimization problem is mainly caused by binary variables. Thus, the method divides the aforementioned operation variables (i.e. decision variables (ii)-(v)) in two categories: (A) operation variables related to binary variables, and (B) other operation variables not related to binary variables. Group A includes the operation variables of those technologies characterized by an on/off status, i.e. the on/off status and the input/output power of heat pump, micro gas turbine, fuel cells and electrolyzer (subset M A ). Since such variables represent the main source of computational complexity, their number is limited by using typical design days. Group B includes all other decision variables, namely the imported/exported grid energy, the charging/discharging energy and the stored energy, and the power generated by those conversion technologies where the on/off status is neglected, i.e. renewable-based technologies and boiler (subset M B ). These variables are defined at every hour of the year.
This classification of the decision variables affect the formulation of both constraints and objective function. In particular, as group A is modeled using design days (e.g. Eqs. (30) and (31)) and group B using every hour of the year (e.g. Eqs. (18)-(22)), the sequence of design days σ couples the two types of decision variables within the energy balances and the objective function (following the same logic of Eq. (36)). The energy balances, Eq. (23), become:  This method features a higher flexibility for the multi-energy system compared to the previous approaches, which comes at the cost of a higher number of decision variables and constraints. Finally, note that for both M1 and M2 σ enables the coupling of the dynamics of conversion technologies of group A during two successive days of the year. The overall optimization framework considered in this work is summarized in Fig. 3. First, the set of typical design days is determined as described above. Based on the adoption of such design days, the optimal design of the multi-energy system is carried out, i.e. technology selection, sizing and operation. Afterward, with the aim of assessing the quality of the different methods, the system design is tested by running the operation with the real hour-by-hour input data, i.e. no design days are used and the system topology is fixed by the design procedure.

Results
In this section, first, the validity of the proposed techniques is tested by considering a simple multi-energy system that can be solved in a reasonable computational time without resorting to design days. Next, the developed tool is applied to minimize the total annual cost and emissions of a multi-energy system characterized by a wider park of conversion and storage technologies.

Comparison of M0, M1 and M2 for a simple MES
The set of available technologies considered in this simplified case is reported on the left-hand side of Fig. 4, which also shows the selected configuration on the right-hand side. The system can include PV panels, boiler, battery and PtG system. The PtG system consists of a PEM electrolyzer (PEME) and a PEM fuel cell using H 2 and air to generate electricity and heat; moreover, it is coupled with a hydrogen tank (HS), which ensures the possibility of seasonal storage. The optimization problem is formulated as a single-objective MILP that minimizes the total annual cost of the system, with no additional constraints on the CO 2 emissions. In this context, a maximum size of 1 MW is imposed for the boiler to ensure the operation of the PtG, which would not be convenient otherwise. Here, the new methods, M1 and M2, are compared with the traditional approach (M0) and the full scale optimization (FSO). The comparison is performed in terms of computational time, system design, operation and annual costs. The MILP is solved by using the commercial software IBM CPLEX 12.7 [61], set to have a relative MIP gap of 0.01% and an aggressive policy for Gomory fractional cuts and probing. Fig. 5 shows the size of the HS and the PEME as function of the design days, which vary from a minimum of 3 to a maximum 72. The PV panels and the battery are not selected, whereas the size of the boiler is constant and equal to 1 MW, i.e. the maximum allowed value. It can be noted that, while the traditional approach M0 is unable to exploit a long-term operation of the HS, both M1 and M2 design a storage size similar to FSO. Indeed, the seasonal operation of the storage calls for a larger size of the HS, as this is charged or discharged during periods longer than one day. The figure shows that, independently of the number of design days, M0 provides a HS size much smaller and a PEME much larger than the one obtained with the FSO. In this case, the size of the storage is set by the daily variations in energy demands and prices, whereas the size of the electrolyzer is set by the maximum value of the thermal demand. On the contrary, both M1 and M2 determine a size of the system similar to FSO when implementing a number of typical days above a certain threshold, which is here of about 25 design days. With respect to M1, a better (less conservative) description of the energy demands is obtained when Fig. 4. Schematic representation of the investigated multi-energy system for model validation: set of available technologies (left) and selected configuration (right). increasing the number of design days, which translates in a more accurate (less conservative) system design. With respect to M2, more design days enable more operation modes of the conversion technologies, which yields a more accurate system operation and design. It is worth noting that M2 is unfeasible when using few typical design days (in this case less than 6), because the produced hydrogen can be reconverted to electricity and heat only according to a limited number of operational modes (i.e. the number of design days), which does not allow to match the thermal supply and demand. A detailed operation of the HS for M0, M1, M2 and FSO during ten days of the year (10.01-20.01) and as obtained in the design phase with 48 design days, is illustrated in Fig. 6. Notably, M0 results in a daily operation policy, which is different for different design days. On the contrary, M1 and M2 result in a long-term operation policy, similarly to the FSO. In the last four days of the examined period shown in the figure, the storage system is operated following an almost-daily pattern. In this case, all methods behave similarly, as highlighted in the insets at the top-right corner of the figure. Overall, the amount of stored energy along the year is equal to 2.7 GWh for M0, 419 GWh for M1, 386 GWh for M2 and 395 GWh for FSO. This confirms a long-term utilization of the storage system by M1 and M2.
To further test the quality of the different methods, full scale optimizations with unfiltered input data and fixed designs, as described in Fig. 3 are carried out. Results are shown in Fig. 7 in terms of total, capital and operational costs as function of the design days. Again, the full scale optimization can be regarded as benchmark as it achieves the minimum objective function. Notably, the trends of the objective function reflect quite well those obtained for the system size: the system design determined by M0, independent of the number of design days, translates into a constant discrepancy in the objective function (around 9% for the total annual cost). On the contrary, a more accurate value of the objective function is determined by M1 and M2 when increasing the number of design days, with the latter approaching faster the value determined by the FSO: a relative error below 1% is obtained by using M2 for a number of design days greater than 24.
Finally, Fig. 8 compares the computation time for the three strategies for a number of typical days ranging from 3 to 72. First, it can be noted that resorting to design days allows to significantly reduce the computation time compared to the FSO, independently of the time horizon description method. Whereas the full scale optimization   requires about 23 h to complete, all the investigated techniques require less than 1 h. Moreover, it can be noted that the additional computation time required by M1 and M2 with respect to M0 is limited. Finally, while M1 has proven to perform generally worse than M2, it however (i) enables the problem resolution with smaller design days and (ii) features a shorter computation time.

Design of a MES of interest
In the following, we aim at designing a MES starting from a comprehensive set of technologies. To this end, we rely on the use of M2, which proved to be the most accurate approximation of the full scale optimization, when based on 48 design days; this proved to be a reasonable trade-off between results reliability and computation time. No additional constraints on the size of the installed technologies are imposed. In this case, the full scale optimization cannot be completed within five days for a MIP gap of 1%, while the optimization problem is solved within 12 h when implementing M2 with 48 design days. The same solver options described in Section 3.1 are used, but for a MIP gap of 1%. The ε-constraint method is used to minimize both total annual cost and CO 2 emissions of the system: first, two single-objective optimization problems minimizing the total annual cost and the total annual emission of the system are solved, providing the upper and lower limits of annual CO 2 emission, respectively. Next, the emission interval is divided into ten steps and single-objective optimization problems are solved to minimize the total annual cost of the system subject to the corresponding emission constraint. Fig. 9 shows the available technologies on the left-hand side and the selected configuration on the righthand side: while the system can consist of several conversion and storage technologies, all design points that lay on the cost-CO 2 emission Pareto considers (i) PV and TS panels, (ii) electrical-driven heat pump, (iii) boiler, (iv) battery, (v) sensible thermal storage and (vi) PtG as optimal technologies. These technologies are selected and sized differently for the different Pareto points. Fig. 10 reports the Pareto front of a reference configuration featuring an available area for solar installation equal to 6000 m 2 , which is based on the actual availability of the investigated neighborhood, and highlights the capital and operational cost contributions. It is worth noting that, on the one hand, a more expensive system design is required when reducing the carbon emissions because of the higher capital cost of solar panels and storage systems. Especially, two sharp slope variations, which are due to changes in the hub topology, can be identified when reducing the CO 2 emissions: at about 400 ton /yr CO2 a significant reduction in the boiler size is observed along with a higher utilization of the thermal storage; at about 280 ton /yr CO2 PtG and batteries are added to thermal storage. Before the former, the decrease in CO 2 emissions does not imply a significant increase in cost; after the latter, the decrease in CO 2 emissions comes at a very large increase in the total annual cost. On the other hand, such designs allow to reduce the operation costs substantially, as less natural gas and electricity are imported from the grids.
Moreover, Fig. 10 compares the Pareto set with a conventional scenario (CS) where electricity is bought from the grid and heat is generated using a boiler. Notably, there are design configurations for which a potential reduction in both cost and emissions can be achieved: taking as an example the configuration at 340 ton /yr CO2 (characterized by a 50% emission reduction with respect to the maximum achievable emission reduction), a cost and emission reduction around 22% and 73%, respectively, are achieved with respect to CS.
To investigate the influence of renewables availability on the Pareto  front, a sensitivity analysis is performed on the area available for solar installation A. Findings are shown in Fig. 11 which presents the Pareto sets for three scenarios characterized by A = 3,000, 6,000 and 12,000 m 2 , i.e. half, reference and double available area, respectively. As expected, a higher renewable installation allows to import less energy from the electrical grid and to achieve lower CO 2 emissions. However, in this case larger storage facilities are needed to fully exploit the higher amount of renewable excess generation, thus translating into higher total annual costs at low emission levels (see Fig. 12). Furthermore, it is possible to notice that the three Pareto fronts coincide at high CO 2 emissions, whereas they are separated at low emissions, where the storage technologies play a relevant role.
A better insight on the different hub topologies is illustrated in Fig. 12, which reports the size of the installed technologies as a function of the CO 2 emission reduction. The CO 2 emission reduction is equal to 0 for the minimum cost configuration and 1 for the minimum emission configuration. A few considerations can be made: i. Natural gas based fuel cells and micro gas turbines, are never installed, independently of the required reduction in CO 2 emissions. ii. The area available for solar installation determines the minimum carbon emissions that can be reached by the multi-energy system. Independently of its value, the total available area is used uniquely for PV panels until about a 80% emission reduction, whereas thermal solar panels are installed for a further emission reduction thanks to their higher conversion efficiency ( Fig. 12a and b). iii. Different combinations of boilers and heat pumps are installed along the Pareto fronts to satisfy the thermal demand ( Fig. 12c and  d), with the former being used only at very low emission reductions (lower than 10-20%). The higher the area available for solar installation, the lower the size of the boiler. Interestingly, a significant heat pump size is always selected despite the high installation cost. This can be explained through an approximate calculation based on average energy prices and conversion efficiencies. Assume the minimum-cost scenario, with no renewable installations. Moreover, assume a cold winter day with conditions favorable to the boiler: average gas price of 0.06 €/kWh, average electricity price of 0.12 €/kWh, average boiler efficiency of 0.92 and average heat pump coefficient of performance of 2.3 (calculated for an ambient temperature of −15°C). Dividing the average energy prices by the average efficiencies leads to an average thermal energy cost of 0.065 €/kWh for the boiler and of 0.052 €/kWh for the HP. A larger discrepancy is observed when repeating the calculation for warmer days of the year. Although approximate, this calculation shows an operational advantage of the heat pump over the boiler, which offsets the higher unit capital cost. Even greater benefits provided by heat pumps are observed in terms of carbon emissions. Consider a specific emission coefficients of 0.137 and 0.237 ton /MWh CO2 for electricity and gas, respectively. Dividing these emission coefficients by the average efficiencies results in an average emission of 0.258 and 0.060 ton /MWh CO2 for the thermal energy generated by the boiler and the heat pump, respectively. Finally, note that the economic and environmental impacts of heat pumps are further reduced when using solar electricity, suggesting a significant potential for the electrification of the energy system under investigation (no natural gas is imported for emission reductions greater than 10-20%). A similar logic explains the absence of fuel cells and micro gas turbines. iv. With respect to the storage technologies ( Fig. 12e-g), thermal storage is installed along the entire Pareto set, with a larger size for lower CO 2 emissions and higher renewable installation. Note that the thermal storage system also helps to satisfy the thermal demand below the heat pump operation limit (commercial heat pumps can typically operate between 10% and 100% of the rated power). On the contrary, battery and power-to-gas systems are only selected at high emission reductions, due to the necessity to fully exploit the renewable generation. Going from low cost to low emissions, the battery is the first to be installed due to the lower installation cost and the higher round-trip efficiency on short-term periods. Then, at very high emission reductions (greater than 80-90%), where a seasonal compensation of the renewable generation is necessary, the PtG system is installed due to higher round-trip efficiency on long-term periods. v Independently of the area available for solar generation, the electrical grid plays a relevant role along the entire Pareto front (Fig. 12h). For the base scenario, even when the maximum solar area is installed (minimum-emission with A = 6,000 m 2 ) the imported electricity still accounts for more than half of the overall energy required. This highlights the potential for reduced grid emissions, which is in fact fairly low within the Swiss framework. In contexts characterized by a higher environmental impact of the electrical grid, e.g. US or UK, the decarbonization of the electrical grid could reduce the need for decentralized storage or the necessity for grid-independent energy hubs. The theme of energy autarky in residential applications has been investigated, among others, by McKenna et al. in [62].
The yearly operation of the storage technologies is further investigated in Fig. 13. In the top sub-figures, the average hourly value of the total stored energy is reported during each month of the year for (i) the minimum-cost design, (ii) the minimum-emission design, (iii) 50% emission reduction and (iv) 90% emission reduction. The total stored energy is given by the sum of the energy stored in the battery, thermal storage and hydrogen storage systems. In the bottom sub-figures, the contributions of the different storage systems are reported with hour resolution for the minimum-emission scenario. It is worth noticing that the seasonal compensation of renewable generation is not required for emission reductions lower than 90%. In fact, a lower amount of energy is stored in summer until a 60% emission reduction due to the lower thermal demand. Thus, shorter operation cycles, though longer than one week, are implemented in this case by means of the thermal storage system. On the contrary, seasonal operation cycles become optimal for emission reductions above 90%. In this case, the optimal storage system to enable seasonal compensation is the PtG due to the low energy losses of the HS, whereas the battery and the thermal storage system are used for shorter compensation, where their higher round-trip efficiency is rewarded. This difference in storage utilization translates into a maximum installed size for the HS significantly bigger than for HWTS and LiB, though none of the storage technologies reaches the maximum size value (see Table 2) independently of the considered scenarios. Moreover, it is observed that doubling the availability of renewable sources a remarkable increase in seasonal compensation is observed, with the HS growing more than three times with respect to the reference case.
Furthermore, sensitivity analysis is performed to investigate the influence of user demands on the design of the integrated system. In particular, while the overall required energy is kept constant, as well as the shape of the demands, three different values of the ratio between the maximum thermal demand and the maximum electrical demand are considered, namely 2.5, 5 (reference value based on input data) and 10. The results of the analysis are reported in Fig. 14, which shows the Pareto sets (left) and the size of the hydrogen storage (right) for the three cases. On the one hand, it is possible to notice that a higher thermal-to-electrical demand ratio leads to lower CO 2 emissions. This is due to the possibility of converting electrical energy into thermal energy with a heat pump that features an average COP of 3.5 (typical value based on average temperature conditions). On the other hand, this emission reduction comes with a higher total annual cost. This depends on the necessity of installing a power to gas system to store the excess in renewable generation. In fact, the power to gas system is only installed above a certain thermal-to-electrical demand ratio due to the different yearly distributions of the two demands (see Appendix): whereas the electrical demand is fairly constant along the year, the thermal demand exhibits a strong seasonal variation. Seasonal operation cycles are required when a high offset between renewable generation (higher in summer) and demand (higher in winter) is registered, i.e. when the thermal demand is predominant over the electrical demand.
Finally, a sensitivity analysis is conducted on the installation cost of battery and power-to-gas to investigate the impact of a cost reduction  on the design of the integrated system. In specific, a cost reduction of 25% and 50% is considered for both technologies with respect to the reference scenario. The cost of the hot water thermal storage is not changed, being this a mature technology. On the one side, reducing the cost of the PtG system does not affect the optimal MES design, as the PtG is always installed at high emission reductions and sized based on the variation between renewable generation and user demand. On the other side, reducing the cost of the battery affects the design of the MES. Results are reported in Fig. 15, showing the size of the battery (left) and the size of the hydrogen storage (right) for the three cases. A reduction in the installation cost of the battery translates into a larger battery installation and a smaller size of the HS, though a certain size of HS is always needed to compensate seasonal variations. Fig. 16 (left) reports the operation of LiB and HS systems along the year for the minimum-emission scenario and a 50% cost reduction. It is observed that although the PtG system represents the higher contribution to seasonal operation, the battery starts being operated with longer-term cycles: from daily cycles during the winter period to weekly and biweekly cycles (up to 400 h) during the spring and summer periods, characterized by a higher variation in renewable generation. This operation of the battery system is justified by a round-trip efficiency significantly higher than that of PtG during bi-weekly periods, as illustrated in Fig. 16 (right).

Conclusions
This paper presents an optimization framework for the optimal design and operation of multi-energy systems including seasonal energy storage. Two novel MILP models are formulated based on the coupling of typical design days and on the distinction of the decision variables between related and unrelated to binary variables. The proposed approaches allow to consider a year time horizon with hour resolution whilst significantly reducing the computation complexity of the optimization problem.
First, the validity of the proposed techniques are tested by considering a simple multi-energy system that can be solved in a reasonable computation time without resorting to typical days, i.e. through a full scale optimization. In this context, a minimum-cost optimization problem is formulated and the proposed strategies are compared against the full scale optimization and the traditional approach relying on uncoupled design days. Findings show that the proposed approaches provide results in good agreement with the full scale optimization, enabling the design, operation and analysis of complex multi-energy systems involving long-term storage.
Furthermore, the developed tool is applied to perform the multiobjective optimization of a residential multi-energy system including seasonal storage. The results are presented in terms of cost-emission Pareto fronts and underline the possibility of a significant reduction in total annual cost and emissions with respect of traditional systems. The topology of the system along the Pareto sets is presented to reveal the behavior of the system and to identify the potential of the investigated technologies for different emission levels and for different boundary conditions. In particular, it is observed that seasonal operation cycles become convenient when a significant reduction in CO 2 emission is required. Moreover, they are favored when a large amount of renewable generation is available and when the system is characterized by a large ratio of thermal to electrical demand. Finally, a remarkable potential for the electrification of the system under investigation, as well as a relevant role of the electricity grid, are observed. the natural gas price, the electrical demand, the thermal demand, are reported in Fig. 17.