Electrochemical conversion technologies for optimal design of decentralized multi-energy systems_ Modeling framework and technology assessment

A R T I C L E I N F O


Introduction
During the 20th century, the energy industry succeeded in the difficult task of meeting the rising electricity and heat demand, which was driven by growing populations and economies. Indeed, this was possible thanks to the use of fossil fuels in large-scale centralized plants, which supply power through national and internationals transmission and distribution grids. However, the evidence that the anthropogenic alteration of the earth carbon balance is leading to climate change has triggered the necessity of finding new routes for energy provision, where no-carbon emission is a must-have feature. In this context, energy efficiency and renewables are the most important building blocks of a pathway consistent with the COP 21 goal of keeping global warming "well below" 2°C [1]. It is well recognized by both institutional bodies and private organizations, that distributed energy technologies represent a key enabler of these measures, as they open significant opportunities to reduce carbon emissions and improve efficiency, while keeping the costs limited, and retaining system reliability [2][3][4][5][6][7]. As a result, both the energy market and infrastructures are being transformed by the development of decentralized energy systems. Especially, the benefits of going decentralized on the energy efficiency and renewable penetration are maximized when multi-energy systems (MES), which exploit the interaction between different energy carriers (e.g. electricity, natural gas and heat), are adopted [8]. On the other hand, the complexity of such systems, which is largely dependent on (i) the number of technologies that can be combined, (ii) the different operation modes that can be adopted, (iii) the uncertainty linked to non-dispatchable generation (i.e. PV and wind) and energy demands, and (iv) the topology of multiple MES integrated with the distribution grid, make the design and operation of such systems particularly challenging. As a consequence, this has prompted many researchers to develop algorithms and computer tools that can tackle such problem. For an overview of the relevant activity in the field, readers are referred to the works by Keirstead et al. [9] and Allegrini et al. [10], who reviewed simulation models within the framework of distributed multienergy systems. In particular, mixed-integer linear programming (MILP) has been particularly favored as optimization method to design and operate multi-energy systems thanks to its flexibility in reproducing complex systems while keeping the computational effort limited. Accordingly, the number of works on this topic is significant. Relevant analyzes focused on the optimal sizing and operation of the energy system (e.g. [11][12][13][14][15][16][17]), on the multi-objective nature of the problem (see e.g. [18][19][20][21]), and on energy networks (see e.g. [13,[22][23][24][25][26]). Along with this, but in a limited number of examples, experimental tests research is carried out (e.g. [27]). However, in the majority of aforementioned works a simplified description of the technology conversion performance is implemented, which neglects the change of the performance with the size and the partial-load operation. Moreover, the dynamics of the conversion technologies is often limited to the start-up/ shut-down or minimum run-time (e.g. [13,28]). The reasons behind this simplified approach lie in the complexity of the optimization problem, but also in the lack of a bridge between the research at the technology/ process level and the research on optimization of MES. This has hindered the adoption of models, more or less detailed, for the conversion technologies, which are in most of the works represented with a constant efficiency in the matrix of coefficients, or, in a few cases, by a set of given efficiencies for part load operation (the description and formulation of the input/output relationships in an energy hub via matrix modeling, with and without constant efficiency, is comprehensively reported in Chicco and Mancarella [29] and Mancarella et al. [30]). Notably, an alternative approach to the constant efficiency formulation was already described by Bloomfield and Fisk in their pioneering work of 1981 [31], where they adopted a piecewise function for the description of heat pumps in a linear optimization problem. More recently, a number of researchers started to tackle this problem within the energy hub framework: Salgado and Pedrero showed how the performance of a combined heat and power (CHP) plant, where the heat and electric output are not independent, can be formulated within an MILP problem [32]; Evins et al. have adopted an improved version of the Bloomfield and Fisk's algorithm for better description of heat pumps operating at partial load [28]; Zhou et al. studied the impact of offdesign performance of boilers and engines on the optimal design and operation of combined cooling, heating and power systems [33]; Bischi et al. developed a piecewise affine approximation for modeling the partial-load performance of the conversion units (boilers, heat pumps and refrigerators) in a MILP [34]; Bracco et al. distinguished between two different values for the maximum and rated efficiencies of some conversion units [35]; Finally, Milan et al. investigated the nonlinear performance of fuel cells, micro-turbines, and internal combustion engines within a MILP [36].
In the framework of multi-energy systems, electrochemical conversion technologies, namely fuel cells and electrolyzers, are regarded as key elements of the technology portfolio. On the one hand, fuel cells (FC) are electrochemical devices that allow for co-generating electricity and heat with very high electrical efficiency regardless of the technology size. Moreover, thanks to their modularity and to the different types of electrolytes that can be adopted, FC can be deployed both at micro and large scale (i.e. kW and MW scale, respectively) without affecting the performance significantly. On the other hand, electrolyzers, which generate hydrogen and oxygen by absorbing electric power through the splitting of deionized water, constitute one of the few technologies that enable seasonal storage of energy with limited (virtually zero) losses of energy over time. Despite the vast research on the use of electrochemical devices in multi-energy systems, there is no clear framework for a reliable description of these technologies in the formulation of such optimization problems. Fuel cells and electrolyzers are either described with constant performance [15,20,37] or as black box models based on manufacturer or literature data [36,38,39]. In all cases, the performance is not clearly connected to first principle thermodynamic models. This is particularly important whenever an integrated system perspective is adopted to provide guidelines for the development of new devices or to compare different fuel cells.
With this contribution we want to overcome these shortcomings, and in particular we aim at (i) providing a clear methodology for the computation within MES optimization problems of realistic electrochemical performance, which are based on first principle models and which can be tuned for whatever electrochemical device; (ii) providing reliable reduced order models for simulating electrochemical devices within the optimization framework, that can be easily implemented in future research works; (iii) assessing the impact of various modeling approximations on the optimal design of integrated multi-energy systems that include electrochemical devices, i.e. determining the most suitable level of detail for simulating these technologies from an integrated system perspective; (iv) making use of these new tools to identify key guidelines for the adoption of different fuel cells in multienergy systems. To this end, first-principle models are implemented to determine the performance of different electrochemical devices under various operating conditions. These are then linearized using different approaches and implemented in the MILP formulation. The optimization framework presented in [40] is adopted to evaluate the impact of approximate conversion efficiency and conversion dynamics on the minimum-cost design of the integrated system. The methodology presented in this work has been applied to the main commercial electrochemical devices, including: proton exchange membrane fuel cells (PEMFC), solid oxide fuel cells (SOFC), molten carbonate fuel cells (MCFC) and polymer membrane electrolyzers (PEME). Note that this methodology can be easily applied to electrochemical technologies not included in this work and holds when considering different objective functions of the optimization problem (e.g. minimum CO 2 emissions). The modeling framework for electrochemical devices presented in details in this work has recently been adopted by Murray et al. [41], and Gabrielli et al. [40].
The paper is structured as follows. Section 2 describes the methodology implemented to model the investigated devices, from the thermodynamic models to their embedding in the optimization framework. Section 3 introduces the optimization framework, whereas results and discussion are presented in Section 4. Finally, in Section 5 conclusions are drawn.

Modeling of conversion technologies
The following electrochemical technologies are considered in this work: • a solid oxide fuel cell (SOFC) using natural gas (NG) as fuel and air as oxidant; • a polymer exchange membrane fuel cell (PEMFC) using NG as fuel and air as oxidant; • a molten carbonate fuel cell (MCFC) using NG as fuel and air as oxidant; • a PEMFC using H 2 as a fuel and air as oxidant; • a high-pressure PEM electrolyzer (PEME) generating hydrogen and oxygen.
These technologies are selected as they represent promising options, already at a commercial level, for decentralized cogeneration and hydrogen production. Other technologies have not been considered because not suitable for decentralized applications or not yet at a commercial level (e.g. solid oxide electrolyzer -SOEC). The last two technologies, i.e. the H 2 -PEMFC and the PEME, are part of a power-togas (PtG) system.
For every technology considered in this work, a systematic modeling procedure is implemented as shown in Fig. 1. First (step 1.a), the polarization curves of the different technologies are reproduced in Matlab® with first-principle lumped models based on physical equations and accounting for electrochemical features and transport phenomena (mass transfer) within the cell [42]. This allows for a proper description of the behavior of the cell stack, typically nonlinear. Second (step 1.b), a wider process simulation is implemented in Aspen Plus®, which not only includes the output of the stack model, but also implements energy and mass balances for all components present in the fuel cell devices (e.g. fuel processor, pumps, compressors, heat exchangers). This allows to predict the on-and off-design performance of the different units. Next (step 1.c), all the considered technologies are simulated in Matlab-Simulink® to investigate their dynamic behavior in terms of generated power (but for the MCFC, whose dynamics is taken from literature). For the PEMFC and the SOFC, this step-methodology builds up and complements the work of Barelli et al. in [43][44][45]. On the contrary, the electrolyzer consists of a simpler layout, and therefore it is simulated without the Aspen Plus® step. A description of the electrochemical and process models developed for this work is reported in Appendix A. In particular, the modeled polarization curves and the corresponding experimental data, as well as the full process layouts for all considered technologies, are reported in Fig. 14, while the main equations implemented to model the electrochemical devices, as well as the main parameters for process simulations, are summarized in Table 5. It is worth noting that the models presented in this work aim at describing commercial products, used in a non-reversible fashion, and optimized to generate electricity and heat (fuel cells) or hydrogen (electrolyzers).
Furthermore, to enable the implementation of these non-linear models within a MILP, the first-principle models need to be reduced to linear models (step 2 in Fig. 1). The linearization approach developed in this work is explained in detail in the next subsection. Finally, the derived reduced order models are implemented within the MILP optimization framework to assess the impact of the aforementioned modeling assumptions and to accurately investigate the performance of the integrated system (step 3 in Fig. 1).
As an exemplary outcome of the thermodynamic models, Fig. 2 shows the performance of the SOFC in terms of electrical power generated by the fuel cell for various partial-load and temperature conditions (left), as well as for a transient of the generated thermal power (right). Such curves are obtained by adopting the simulation procedure described in Appendix A for various operating temperatures, namely 650°C, 750°C (corresponding to the polarization curve in Fig. 14a) 850°C.

Reduced order models -On/off-design performance
The developed first-principle models describe the conversion performance of the considered systems for various operating conditions, e.g. partial-load, temperature, pressure. At design phase, we assume that reference operating conditions are maintained at any power level through the device control system. Therefore, only the partial-load dependency of the conversion performance is considered in the following. This dependency, typically nonlinear (see Fig. 2), is linearized for use in an optimization framework. In particular, for each of the considered conversion technologies the outlet power P is expressed as an affine function of the inlet power F and of the technology size S, i.e. of the rated input power. For the fuel cells, P and F indicate the generated electrical power and the inlet fuel power (fuel LHV), respectively. For the electrolyzer, P and F indicate the generated hydrogen power (hydrogen LHV) and the inlet electrical power, respectively.
As for the linear models, we consider piecewise affine (PWA), affine (A1), and linear (A2) approximations. Such approximations and the corresponding errors with respect to the detailed thermodynamic models are illustrated in Fig. 3 for the SOFC, NG-PEMFC, MCFC and PEME. The affine approximation A1 is the convex hull of the original curve, whereas the linear approximation A2 is the best fit line which goes through the origin. The PWA approximation implements different affine segments for different sections of the performance curves. The optimal position of n breakpoints is determined through the nonlinear optimization problem described in [42]. Differently from [42], the first and last breakpoints must coincide with the first and last point of the original curve, respectively. Here, 5 breakpoints, i.e. 4 line segments, are used for all the considered technologies, as this proves to describe well the exact performance.
The following equations express the output power for the A1 (Eq. (1)), A2 (Eq. (2)) and PWA (Eq. (3)) approximations: where for the A1 and PWA approximations while for the A2 approximation Here, the parameters α β γ , , express the power and size dependencies of the generated electrical power, with the index i referring to the i-th affine segment of the PWA approximation; the parameters δ ζ κ ν , , , are used to constrain the input power between its minimum and maximum values; S is the size of the technology; ∈ x {0,1} is a binary variable specifying whether the device is turned on, thus producing power but also incurring in the auxiliary consumption + δS ζ . It is worth noting that due to the concavity of the curve, no additional integer variables are required in the PWA approximation to select the active line segment.
Fuel cells are cogenerative systems, i.e. they generate electricity and heat at the same time. Therefore, linear correlations are developed for the generated heat, likewise for the electrical power. As the generated thermal power Q shows a small deviation from linearity, only the affine approximation A1 and the linear approximation A2 are considered. Such approximations and the corresponding errors with respect to the thermodynamic models are illustrated in Fig. 4 for the SOFC, NG-PEMFC and MCFC. In the two cases, the thermal power is expressed, respectively, as where the parameters α β γ , , T T T characterize the power and size dependencies of the generated thermal power. Furthermore, it must be mentioned that due to the modular nature of the investigated electrochemical technologies, their performance is not significantly affected by their size, thus translating into = = = γ ζ ν 0. All the parameters involved in the description of the conversion performance are reported in Table 1.

Reduced order models -Conversion dynamics
The dynamic behavior of the aforementioned technologies is considered in terms of (i) start-up/shut-down power trajectories, (ii) rampup/ramp-down limitations, and (iii) generation transient behavior. Start-up and shut-down limitations are modeled based on the set of linear constraints proposed by Arroyo and Conejo in [46] (in particular, we refer to Eqs. (1)- (12) in that paper). In addition, a time constant σ based on manufacturer and literature data is introduced in this work to handle start-up and shut-down operations shorter than the length of the considered time intervals, i.e. one hour here. Thus, the input powers at the h-th period of the start-up and shut-down processes are constrained, respectively, as where τ SU and τ SD indicate the duration of the start-up and shut-down processes, respectively; R RU and R RD indicate the ramp-up and rampdown limits, respectively; F min indicates the minimum power. Furthermore, the number of start-ups in one year is limited by a maximum value N SU : t is a binary variable indicating whether the unit is started-up during the t-th time instant, with T the length of the time horizon.
The generation transients are modeled through a first-order dynamics that approximate the nonlinear behavior of the considered technologies. As an example, Fig. 5 shows the nonlinear dynamics of the SOFC when varying the amount of thermal power provided to the end users, as well as the corresponding first-order approximation. Such a linear approximation, i.e. the red line in the figure, is determined by matching the overall amount of energy provided by the thermodynamic model, i.e. the blue and red areas underlying the original and approximate curves. When varying the generated power from q 1 to q 2 , the thermal power follows this first-order, non-homogeneous dynamics: Here, τ is the time constant characterizing the power transient. The relative error between the energy provided by the thermodynamic model and that provided by the first-order approximation is about 0.3% over one hour (time interval implemented for the optimization described in Section 3). However, for the considered technologies the time constant τ is typically shorter than the an hour. Therefore, based on the linear approximation a parameter ω is defined that expresses the fraction of energy lost or gained during the power transient between two consecutive hourly time instants. In this context, the power is assumed to change with an hourly basis and the parameter ω indicates the amount of energy lost or gained during an hour. Assuming = ω 0 corresponds to neglecting the power transients, i.e. going instantaneously from one power level to another. Thus, for all time steps ∈ … t T {1, , } the generated energy E t is calculated as All the parameters involved in the description of the system dynamics are reported in Table 2.

Formulation of the optimization framework
To evaluate the different methodologies introduced in Section 2, the various technologies are considered within the framework of the multienergy system (MES) presented in Fig. 6. The considered MES is a set of conversion and storage technologies that has the primary goal of supplying electrical and thermal energy to a residential neighborhood in Zurich, Switzerland. The MES is connected to the natural gas and electrical grids, and the energy can be converted locally through all the technologies described above, namely (i) SOFC, (ii) NG-PEMFC, (iii) MCFC, (iv) H 2 -PEMFC, (v) PEME. The generation and re-conversion of hydrogen through the PEME and H 2 -PEMFC is referred to as power-togas (PtG). The two devices are connected to a hydrogen tank (HS) which stores hydrogen and decouples its production and utilization. Furthermore, a hot water sensible thermal storage (HWTS) is utilized to satisfy the thermal demand when this decreases below the minimum value achievable with the fuel cells. In this work, an electrical storage is not considered because beyond the scope of this paper. Readers can refer to Gabrielli et al. [40] for a comparison among different storage systems via the present modeling framework. Although gas and electricity can be imported/exported, no network models are considered in this work, as the focus is on the local energy conversion and storage within the boundary of the considered MES.
A MILP is formulated to determine the selection, size and operation over one hour. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.) of the considered conversion and storage technologies, based on the minimum cost of the integrated system. In particular, the input data of the optimization problem are ambient temperature, the prices of electricity and gas, the electrical and thermal demand profiles, and the set of available conversion technologies with the corresponding performance and cost coefficients. Outputs of the optimization problem are the selection and size of the installed technologies, the scheduling (on/ off status) of the conversion technologies, the input and output power of the conversion technologies, the stored, charged and discharged energy of the storage technologies, and the imported/exported power from/to the electricity and gas grids. The MILP can be written in general form as where c and d represent the cost vectors associated to continuous and binary decision variables, v and w, respectively; A and B are the corresponding constraint matrices and b is the constraint known-term; N v and N w indicate the dimension of v and w, respectively. The binary variables are introduced to model the operation (on/off) and the performance and costs of the considered technologies. Within this framework, both continuous and binary variables are optimized, with the latter being introduced to model performance and costs of the investigated technologies. In the following, the optimization problem is described in terms of input data, decision variables, constraints, and objective function. A detailed description of these features can be also found in Ref. [40].

Input data
The input data are time-dependent profiles for 2016 with an hour resolution in Altstetten, a neighborhood of Zurich. The neighborhood is characterized by a peak electricity and thermal demands of 0.43 and 2 MW, respectively, and a total annual electricity and heat demands of about 1.6 and 4 GW h, respectively. The hourly demand data are determined and discussed by Murray et al. [41], and are reported in Appendix A of Ref. [40]. Note that while clusters of typical design days are generally implemented to model the input data, the detailed hourly profiles are considered here by using the clustering methodology M2 introduced in Ref. [40] for modeling the system operation along the time horizon. These profiles are assumed to be constant along the lifetime of the system and known without uncertainty. In other words, we assume that the possible realizations of the uncertainty are represented by the historical data. Inputs to the optimization problem are: i. The ambient temperature ∈ A T , where T indicates the length of the time horizon. Data are available at every hour of the year, thus = T 8760.
ii. The electricity and heat demands, L e and ∈ L T h , respectively. iii. The electricity and gas prices, u e and ∈ u T g , respectively. While the import gas price is considered to be fixed along the year at 0.064 €/kW h [47], the Swiss electricity spot market price for 2016 (varying from 0.01 to 0.12 €/kW h) is considered [48]. Note that the same price is considered for electricity import and export. iv. The set of available conversion and storage technologies described above with the corresponding performance (Table 1) and cost ( Table 2 in [40]) coefficients.

Decision variables
The following decision variables are obtained as outputs of the optimization problem: i. The size of the installed technologies, ∈ S M , where M indicates the number of the available technologies. Note that determining S (greater than zero if a technology is selected) also implies selecting the technologies. ii. The on/off status of the conversion technologies ∈ The operation of fuel cells and electrolyzer (decision variables ii and iii) is modeled through 48 typical design days, whereas the operation of the storage devices and the imported/exported energy (decision variable iv and v) are determined at every hour of the year. This is done by implementing the M2 method presented in [40].

Constraints
The optimization constraints can be divided into two categories: i. Performance of conversion and storage technologies. The performance of the conversion technologies are described by the reduced order models introduced in Section 2, i.e. by Eqs. (1)- (12). The behavior of the storage technologies is described by the following linear dynamics: Here, Λ and Π are self-discharge parameters; g(A t ) expresses the influence of the ambient temperature on the losses of the storage devices, as suggested in [49] for the thermal storage; η indicates the charging and discharging efficiency, considered to be equal; t Δ is the duration of the time interval t τ ; is the time required to fully charge or discharge the storage. Note that there is no need for binary variables specifying whether the device is charging or discharging. Indeed, if η is lower than one, the optimizer automatically selects one of the two modes due to the losses in the charging/discharging process (according to the MIP gap of the optimization). The periodicity constraint, Eq. (15), imposes the same storage level at the beginning and at the end of the year. The performance coefficients of the storage technologies are reported in Table 1 in [40]. ii. MES energy balances. The considered energy carriers are electricity, heat, natural gas and hydrogen. In particular, natural gas and hydrogen can be used as fuels by the fuel cells to generate electricity and heat; electricity can be used by the electrolyzer to generate hydrogen. The sum of imported and generated power must equal the sum of exported and used power for all energy carriers j, for all time intervals ∈ t T {1, }: Here, i indicates the i-th technology, U 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. Limitations on the imported and exported energy could easily be included when necessary.

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. The equivalent annual investment cost is computed through the annuity factor a, where an interest rate of 6% is considered. The annual operation cost is calculated based on the amount of imported and exported electricity and gas during the year: The annual maintenance cost is given as a fraction ψ of the annual capital cost [50,51]: The cost coefficients for all the considered technologies but the MCFC can be found in Table 2 in [40]. The cost of the MCFC varies between 3000-5000 €/kW of generated electricity and can be found in [51].

Technology modeling definition
The optimization framework is used to investigate the influence of the modeling methodologies on the optimal design of the integrated system. To this end, the following technology modeling (TM) scenarios are proposed: TM I. Reference. The conversion performance is modeled by using the piecewise affine (PWA) approximation. The system dynamics is modeled using the set of linear constraints described in Section 2.2. This methodology represents the most detailed description of the conversion technologies that can be implemented within a linear optimization framework, and it is therefore regarded as the reference scenario. TM II. Linear conversion. The conversion performance is modeled by using the affine approximation A1. The system dynamics is modeled using the set of linear constraints described in Section 2.2. In this way, the effect of a linear description of the conversion performance is evaluated. TM III. No dynamics. The conversion performance is modeled by using the piecewise affine (PWA) approximation. No dynamic behavior is considered. This implies that no start-up/shut-down and ramp-up/ramp-down limitations are accounted for, i.e. the conversion devices can immediately (i) be switched on/off, and (ii) move along the off-design curve instantaneously. This allows to evaluate the importance of considering the dynamic behavior of the conversion technologies. TM IV. Simplest description. The conversion performance is modeled by using the linear approximation A2. No dynamic behavior is considered. This is the simplest possible description of the underlying conversion units, i.e. constant efficiency. No binary variables are required in this case, as both the system dynamics and the minimum-power limitations are neglected.

Evaluation of the technology models
The aforementioned modeling methodologies are evaluated through the following procedure. First, the optimal design is performed for all methods TM I-IV. As TM I is without any doubts the closest to the thermodynamic models, the designs obtained with TM II-TM IV are tested by operating them with the detailed technology description given by TM I (see Section 2). When operated on TM I, designs given by TM II-IV translate into a violation of the energy balances (Eq. (18)). This means that the energy demands can either be unsatisfied or exceeded. In the former case, the energy demand required by the end-user is not met. In the latter case, too much energy is produced (e.g. due to the impossibility of switching off the unit when the storage is already full) and needs to be dissipated to the environment. In specific, while the electricity demand can always be satisfied by directly importing electricity from the grid, the thermal demand can only be provided by the conversion units. Thus, to evaluate the violation/dissipation within the optimization framework, the heat balance must be rewritten, for all t, as where Ω v and Ω d indicate the total thermal energy violated and dissipated, respectively, and are non-negative. Furthermore, a fictitious arbitrarily high cost N is associated to any violation/dissipation event to prevent the optimizer identifying this as optimal solution. Thus, the total annual cost is written as where N is here selected equal to 10 5 €/kW h. In the following, the energy violation and dissipation are expressed as fractions of the annual thermal demand:

Results and discussion
First, we investigate the required level of details for simulating each of the considered technologies within integrated multi-energy systems. To this end, the different methods introduced in Section 3 are compared in terms of trade-off between accuracy and computation time. Next, we provide guidelines for the deployment of fuel cell cogeneration units within multi-energy systems. To do so, the most relevant techno-economic parameters from an integrated system perspective are defined and assessed.

Required level of detail for technology simulation
In this section, the impact of the modeling approximation is quantified for each conversion technology. To this end, the optimization problem is solved by considering a MES consisting of one unit at a time coupled with the thermal storage. The different modeling scenarios TM II-IV are compared against the reference scenario TM I in terms of system design, i.e. size of the installed technologies, total annual cost, energy violation and computation time. The results are summarized in Table 3 for each of the considered technologies. All optimizations are solved by using the commercial software IBM CPLEX 12.7 [52], set up to have a relative MIP gap of 1%. All simulations are run by using 6 cores of a computer featuring a Processor Intel(R) Core(TM) i7-4790 CPU, 3.60 GHz and 8 GB of installed RAM.
SOFC-based MES. Linear conversion performance, i.e. TM II, results in approximately the same FC and HWTS sizes determined by the reference approach TM I. In this case, the size of the system is constrained by the slow SOFC start-up/shut-down dynamics. This forces the SOFC to be constantly switched on and results in the most cost-effective tradeoff between the size of the FC and that of the HWTS to satisfy the thermal loads when they are lower than the minimum generated power. Indeed, when this happens the excess power must be stored in the thermal storage and used afterward. Because of the system sizes, TM I and TM II are characterized by the same capital and operation costs (variations smaller than 1% are considered negligible, being below the MIP gap), where the operation cost is determined by operating the TM II on TM I. In this case, a negligible energy violation is observed during three hours of the year, where the peak thermal demand occurs in winter, whereas no energy dissipation occurs. Fig. 7 shows the energy violation/dissipation across the entire year when operating the TM II-IV SOFC (left-hand side) and PEMFC (right-hand side) designs. Here, the TM II energy violation is indicated by the blue shaded area. A little variation is observed in terms of computation performance. When neglecting the system dynamics, i.e. in the TM III case, a greater size of the FC and a smaller size of the thermal storage are obtained. In this case, the fuel cell is characterized by a greater (fictitious) flexibility, being able to immediately switch on/off. This reduces the necessity for thermal storage compared to the realistic description. In this case, the higher installation cost (+1.3%) is compensated by a lower operation cost (−1.9%), which depends on the possibility of dissipating heat during the summer months (733 h). In terms of computation time, a remarkable reduction is observed due to a reduction in the number of binary variables.
The simplest technology description, i.e. TM IV, results in a slightly smaller size of the FC and a significantly smaller size of the HWTS (−54%). This is mostly due to a high fictitious flexibility of the FC, which can now cover the entire range of thermal demand, further reducing the necessity for the HWTS. Moreover, the HWTS is only used to cover the demand peaks (most cost-effective trade-off), translating into a smaller FC size. This design is characterized by lower capital and installation costs, due to the possibility of not satisfying the thermal demand during the winter peak (21 h) and dissipating energy during the summer months (834 h). Note that this violation represents about 10% of the thermal demand required during the relevant period, which could results into a discomfort for the end-users. On the other hand, the computation effort is drastically reduced (−90.7%) due to the absence of binary variables.
Overall, if a violation/dissipation in the order of 1% is required, the system dynamics can be neglected, whereas the minimum-power constraints should be taken into account, making TM III a suitable option for designing SOFC-based systems. If no violation/dissipation is accepted, TM I or TM II must be adopted.
PEMFC-based MES. As for the SOFC, TM I and TM II result in about the same size. This translates into a negligible discrepancy (lower than 1%) of the total annual cost and in no violation/dissipation when operating TM II with the TM I technology description. On the other hand, it is worth noting that this approximation does not simplify the solution of the optimization problem, but results into a larger computation effort with respect to the base scenario.
The size of the FC is not particularly affected when neglecting the system dynamics, being the PEMFC characterized by faster start-up/ shut-down trajectories. However, a smaller size of the HWTS is required because the FC can switch from one level of generated power to another immediately. In this case, some energy violation occurs during the winter peaks, as shown in Fig. 7. Also, a drastic reduction in computation time is obtained (−89%).
The simplest description scenario results in a significant undersizing of the system (−6.3% and −49.2% for the fuel cell and thermal storage, respectively) due to the absence of minimum-power constraints. This translate into a lower total annual cost (−2.9%), due to the possibility of violating the energy demand during the winter peak (41 h), whereas no energy dissipation occurs (thanks to the fast dynamics of the device). Notably, the energy violation (see Fig. 7) represents the 21% of the thermal demand required during the relevant period, which could results into a discomfort for the end-users. Obviously, a drastic reduction in the computation effort is observed (−99%).
Overall, a simple description seems suitable for designing PEMFCbased multi-energy systems. However, if no violation/dissipation is accepted, TM I or TM II should be adopted.
MCFC-based MES. Considerations similar to those made for the SOFC-based system hold in this case. However, the MCFC features a significantly less pronounced nonlinearity and greater minimum power requirement. This implies a wider feasibility region of the optimization problem when neglecting the system dynamics and therefore a greater computation effort. Thus, if a violation/dissipation lower than 1% is required, TM I or TM II are suitable options for designing MCFC-based systems.
PtG-based MES. Similarly to the PEMFC-based system, assuming linear conversion performance (TM II) and neglecting system dynamics (TM IV) result in approximately the same sizes with respect to TM I. This translates into similar costs and negligible violation/dissipation. On the contrary, the TM IV design features sizes significantly lower than the reference description. Together with the limited flexibility of the PtG (where hydrogen cannot be bought from the grid, in contrast to natural gas, but must be locally generated), this results into energy violation during the peak of thermal demand in winter. It is worth   noting that this violation represents about 20% of the thermal demand required during the period where the violation occurs. Overall, if a violation/dissipation in the order of 1% is required, TM III represents the most suitable modeling option for MES implementing PtG. If no violation/dissipation is accepted, TM I must be adopted for simulating the conversion devices.
Summary. The analysis above assesses the most suitable level of details for simulating the considered technologies within integrated multi-energy systems. Considerations are made on the accuracy and computation effort of the system design. Findings can be summarized as follows: • A linear approximation of the conversion performance (TM II) does not significantly impact the accuracy of the results, but does not translate into significant computation improvement.
• Neglecting the system dynamics (TM III) does not significantly impact the accuracy of the results, and it typically improves the computation performance (but for the MCFC).
• The simplest technology description (TM IV), with constant efficiency, no system dynamics and no minimum-power constraints, impacts the accuracy of the results as it translates into substantial energy violation/dissipation. Since it also results into a drastic reduction of the computation complexity, it should be implemented when little accuracy and fast performance are required.
As an example, when an energy violation/dissipation in the order of 1% is required, TM III should be used for SOFCs and PtG systems, TM II should be used for MCFCs and TM IV should be used for PEMFCs. However, if no violation/dissipation is accepted, TM I must be adopted for all the conversion devices (or TM II for SOFC and MCFC).

Design guidelines for FC implementation
In this section, the optimization problem is solved by considering a MES where all the investigated technologies can be installed. Findings indicate that, for the considered input data used as boundary conditions, the PEMFC is always preferred as a generation unit over the SOFC, the MCFC and the PtG, i.e. PEMFC is the only selected technology together with the thermal storage. On the one hand, the power to gas system is expected to be less convenient than the fuel cells as simple cogenerative units. In fact, no renewable generation is considered here, thus eliminating the need of hydrogen storage [40]. On the other hand, further analysis is required to better understand why and when PEMFCs are preferred over SOFCs or MCFCs, to provide guidelines for the optimal design of these systems within the framework of decentralized generation.
First, a very simplified optimization problem, which can be solved analytically, is formulated to compare the different FCs at a conceptual level. As an example, the comparison between SOFC and PEMFC is reported below. For both conversion technologies the following is assumed: (i) constant conversion efficiencies, η, and no system dynamics (equivalent to the simplest description scenario, TM IV); (ii) constant unit installation costs, c; (iii) same lifetime and durability features for both devices; (iv) thermal demand higher than the electrical demand at every time instant, thus determining the amount of generated power. In this case, the total annual cost of the fuel cells can be written as where the subscript ∈ t T {1, } indicates the t-th instant of the time horizon; η T and η E indicate the thermal and electrical efficiency, respectively; L h and L e indicate the heat and electricity demands, respectively; u g and u e are the natural gas and electricity price, respectively. Here, the first term on the right-hand side is the capital cost, with the size L η / h max T being defined by the maximum thermal demand. The second and third terms on the right-hand side are the operation cost of purchasing gas and selling electricity, respectively, both being driven by the thermal demand. Under these assumptions, the cost of the SOFC is lower than that of the PEMFC when the following condition is fulfilled: where the subscripts S and P refer to SOFC and PEMFC, respectively, and Since for typical cost and efficiency values ≃ − f f 1 2 , it is possible to simplify Eq. (26) as In other words, the SOFC is convenient over the PEMFC when the ratio of unit capital costs Γ is lower than the ratio of thermal efficiencies H, i.e. when the capital cost per unit produced thermal power is lower. This also implies that, assuming a constant first-principle efficiency, the SOFC becomes convenient over the PEMFC at a certain value of the thermal-to-electrical efficiency ratio. Although derived for a simple optimization framework, this condition highlights the importance of thermal efficiency within the framework of distributed generation: here the thermal demand often represents the limiting factor for cogeneration and a thermal efficiency larger than the electrical one is thus favored. The same approach can be applied to the MCFC, where the subscript M replaces S.
Following this simplified, yet straightforward analysis, the trade-off between electrical and thermal generation is investigated in detail through a sensitivity analysis on the input data of the full optimization problem. The sensitivity analysis is summarized in Table 4. A number of parameters are varied for SOFC, MCFC and MES to identify the optimal type of fuel cell for a given multi-energy system. First, the unit installation costs of SOFC and MCFC are varied while maintaining the other parameters unchanged.
The case of the SOFC for scenario TM I is illustrated in Fig. 8, which shows the installed size for SOFC and PEMFC, as well as the total annual cost of the system J as a function of Γ , the ratio of unit installation costs. Note that the unit capital cost of the FCs is a function of the size (see Table 2 in [40]). The unit installation cost of the PEMFC is maintained constant, while that of the SOFC is reduced from 100% to 1% of its current value (moving leftward along the horizontal axis). One can notice that the SOFC starts being installed only when its cost is about 55% of that of the PEMFC, i.e. a threshold cost ratio = * Γ 0.55 is identified below which the SOFC is convenient with respect to the PEMFC. This is mostly due to the lower SOFC thermal efficiency, which requires a size larger than the PEMFC to satisfy the thermal demand.
Next, we have computed the cost threshold for different values of the thermal efficiency while keeping the other parameters unchanged. Findings are illustrated in Fig. 9, which shows the cost threshold * Γ as a function of the thermal efficiency ratio H for scenario TM I (solid line in the figure). Note that the design value of the thermal efficiency is taken for each FC, as this varies with the system operation. The SOFC-PEMFC comparison is reported on the left-hand side, whereas the MCFC-PEMFC comparison is reported on the right-hand side. The thermal efficiency of the PEMFC is maintained constant, while the thermal efficiency of SOFC and MCFC is increased from their current value to the PEMFC value (moving rightward along the horizontal axis). This corresponds to increase the current thermal-to-electrical efficiency ratio in the ranges 0.47-1.14 and 0.44-1.58 for SOFC and MCFC, respectively. Note that a higher thermal efficiency results in a lower electrical efficiency, as the first-principle efficiency is kept constant at 0.85 and 0.74 for SOFC and MCFC, respectively (0.9 for the PEMFC). Two regions are identified: (i) the region below the cost threshold (blue-shaded line), where the SOFC and MCFC are selected, being more convenient than the PEMFC; (ii) the region above the cost threshold (red-shaded area) where the PEMFC is selected, being more convenient than SOFC and MCFC. It is worth noting that, independently of the thermal efficiency, the MCFC is selected only when costing less than the PEMFC. This is due to the combination of slow dynamics and high minimum-power constraints. On the contrary, above a thermal efficiency ratio of 98% (i.e. a thermalto-electrical efficiency ratio of 1.1), the SOFC is selected even if more expensive than the PEMFC, i.e. > * Γ 1. This depends on the lower minimum-power constraints and on the longer lifetime (15 years for the SOFC versus 10 years for MCFC and PEMFC), which offset the slower dynamics. Indeed, the lifetime of the technologies plays a role as a longer lifetime translates into a lower annualized cost. However, it is worth mentioning that, for both SOFC and MCFC the current cost is significantly higher than that of PEMFC (about 150% and 200%, respectively), thus suggesting that a remarkable cost reduction is needed for them to become convenient.
Moreover, the importance of SOFC and MCFC dynamics is investigated by studying the cost-efficiency function when (i) reducing the start-up/shut-down times down to 2 h, (ii) reducing the ramp-up/ ramp-down times down to one hour, and (iii) increasing the number of yearly on/off cycles up to 100. The results are shown in Fig. 10  It is worth mentioning that no significant advantage is observed when varying SU/SD, RU/RD, N SU separately. Indeed, each of these factors would prevent the fuel cells to follow the variations of the thermal load. Therefore, a faster conversion dynamics is required in terms of both start-up/shut-down and ramp-up/ramp-down and should be coupled with a greater number of the thermal cycles per year. Furthermore, a sensitivity analysis is performed on the cost of the thermal storage, which can range from the current value to 1% of the current value. However, no significant impact is observed in this case due to the already low cost of the thermal storage.
Moreover, the cost-efficiency function is investigated for different features of the multi-energy system. First, Fig. 11 shows the cost threshold * Γ as a function of the thermal efficiency ratio H for different thermal and electrical demands required by the end user. While the overall required energy is kept constant, as well as the shape of both demands, the ratio of thermal-to-electrical maximum demand = λ L L / h max e max ranges from 2.5 to 10, with 5 being the reference value.
It is possible to notice that SOFC and MCFC become more convenient with respect to the PEMFC when increasing the value of λ, i.e. for endusers with a higher share of electrical demand, due to a more effective utilization of the generated electrical power.
To get an insight of how these systems would behave in terms of self-sufficiency, the behavior of the cost-efficiency function is analyzed when imposing a maximum limit on the amount of electricity imported from the grid I. Fig. 12 shows the cost threshold * Γ as a function of the thermal efficiency ratio H when a maximum threshold of 100%, 50% and 25% of the reference value is imposed on the amount of imported electricity. One can notice that SOFC and MCFC become more convenient with respect to the PEMFC when constraining the amount of imported electricity. Indeed, a higher electrical efficiency is rewarded in this case. Moreover, a higher relevance of the on-site electrical generation requires all the devices to be switched on more often, thus reducing the impact of the slow system dynamics. Finally, the cost-efficiency function is studied for different electricity prices. Indeed, decentralized MES can be characterized by different tariffs depending on the location, the size of the system, the policy framework, and so forth. Fig. 13 shows the cost threshold * Γ as a function of the thermal efficiency ratio H for three different electricity prices, namely the 2016 Swiss spot market price (reference case), and two constant electricity prices equal to 0.1 and 0.2 €/kW h (meaningful values for residential generation). Note that whereas in the former case electricity can be sold to the grid also at the spot market price, a constant selling price of 0.1 €/kW h is considered for the other two cases. One can note that, for the considered constant electricity prices, the SOFC represents the most convenient solution within the whole range of investigated thermal efficiency ratio, due to the longer lifetime and lower minimum power. In this case, an increase in the thermal efficiency plays a less relevant role and even leads to more expensive solutions for high electricity prices.

Conclusions
This paper defines the required level of detail for simulating electrochemical conversion devices within the optimal design of  decentralized multi-energy systems. Furthermore, it provides design guidelines for the optimal development and deployment of fuel cell cogeneration units within decentralized MES. First, detailed thermodynamic models based on a first-principle approach are implemented to accurately represent the part-load performance and the dynamic features of different types of fuel cells and electrolyzers under various operating conditions. Then, as such nonlinear models are intractable for use in the MILP-based optimization of integrated systems, different linear models are developed. In particular, linear and piecewise linear approximations are implemented to describe the conversion performance, and a first-order approximation is used to describe the system dynamics.
Starting from these reduced order models, several modeling simplifications are investigated to quantify the impact of nonlinear conversion performance, conversion dynamics and minimum-power constraints on the design of integrated systems. More specifically, the proposed considered modeling simplifications are compared against the reference reduced order models to evaluate the trade-off between the accuracy and the computation effort of the optimization problem, thus allowing to define the required level of detail to model electrochemical devices within the framework of decentralized multi-energy systems. Findings show that (i) a linear approximation of the conversion performance (TM II) does not significantly impact the costs and feasibility of the integrated system, but does not simplify the optimization problem with respect to a more accurate PWA approximation, which is therefore generally preferred; (ii) neglecting the system dynamics (TM III) does not significantly affect the system costs and feasibility, and it typically improves the computation performance (but for the MCFC); (iii) the simplest technology description (TM IV), with constant efficiency, no system dynamics and no minimum-power constraints, impacts the accuracy of the results as it translates into energy violation/ dissipation. Since it also results into a drastic reduction of the computation complexity, it should be implemented when little accuracy and fast performance are required.
Furthermore, the optimization framework is exploited to provide guidelines for the deployment of FC cogeneration units within decentralized multi-energy systems. This is done by defining and assessing the most relevant parameters from the perspective of the integrated system design. The technology cost, thermal efficiency and conversion dynamics are found to be the most important techno-economic parameters, and different design regions are identified as functions of such parameters. Overall, findings show that, unless a significant cost reduction for SOFC and MCFC occurs, the PEMFC is most suited for decentralized cogeneration due to its faster dynamics and its higher thermal efficiency. However, when the end-users require a higher fraction of electrical energy or when on-site electrical generation is requested, SOFC and MCFC can be preferred. Similarly, the SOFC represents the most convenient option in residential contexts featuring a constant electricity price. Finally, it is worth noting that the proposed analysis holds when considering minimum CO 2 emissions instead of minimum cost as the objective function of the optimization problem.  Table 5 Key equations and parameters of electrochemical models and process simulations for the considered technologies. ] 180 [53] 180/500 [54] 80 [55]   temperature, and (iii) a preferential oxidation (PrOX) reactor for CO abatement. Additionally, heat exchangers for heating up natural gas and generating steam for reforming and WGS are required; the heat sources consist in the products of the anode-cathode outlets combustion and in the WGS outlets (as for the SOFC, minor amount of natural gas may be added to the combustor under certain operating conditions). In addition to the fuel processor, the PEMFC plant includes a detailed system for humidifying the cathode air to the level required by a safe operation of the polymeric membrane [72]. Finally, the plant includes the air blower, pumps, the cogeneration heat exchanger, and power electronics. The process data are taken from Refs. [73,43]. When the PEMFC is fed directly with H 2 , i.e. when it is part of a power-to-gas system, the process model is very much simplified (no need for fuel processor) and based on the commercial system by Swiss hydrogen (10 kW) [74]. Overall, the dynamic behavior of the process components, e.g. heat exchangers, is modeled following the approach presented by Barelli et al. in [43,44]. It is worth noting that, also in the case of PtG, air is used as oxidant 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 : O 2 2 ratio of 2:1 and consumed by the PEM fuel cell with a H : O 2 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 costs and flow rate limitations. A comparison between these two alternatives is discussed by Büchi et al. in [54].
MCFC process layout: the process layout of the MCFC (Fig. 14(f)) reproduces the FuelCell Energy system (300 kW) [75]. The model of the cell stack is based on the work by Audasso et al. [55], which only accounts for steady-state operation. The polarization curve considers voltage losses across the anode, the electrolyte and the cathode, as reported in Table 5. The open-circuit voltage is evaluated through the Nernst equation, where average partial pressures between inlet and outlet are used; the voltage losses are expressed through empirical correlations. Different CO 2 molar fractions at cathode inlet ranging from 2 to 6% are considered. The values of the parameters describing these contributions are based on the references reported in Table 5, and the resulting polarization curves are validated against the experimental data reported in Ref. [55]. As for the SOFC, the MCFC considered in this work features an internal reforming, thus reducing the fuel processor to an adiabatic pre-reformer and a heat exchanger for steam generation and anode inlet heating. The anode inlet is at 600°C and features the following composition (vol): H 2 0.03, CO 0.01, CH 4 0.27, H O 2 0.69. Notably, the MCFC process layout features a combustor on the cathode inlet line: by burning the anode outlet stream with the cathode outlet stream and fresh air, the partial pressure of CO 2 in the cathode inlet is increased to the desired level for enabling the carbonate ion exchange in the membrane. Finally, the process layout includes pumps, a cogenerative heat exchanger, an air blower, and power electronics. The process data are taken from Ref. [76]. The dynamic features of the device are based on literature data [77,78]. PEME process layout: the process layout of the electrolyzer (Fig. 14(h)), which is based on the commercial system Silyzer by Siemens (100 kW) [79], is relatively simple compared to that of the fuel cells as no fuel conversion is needed. The plant consists of (i) water pumps, (ii) H 2 and O 2 purification systems, (iii) cogenerative heat exchanger, and (iv) power electronics equipment (because of its simplicity, the PEME process was directly modeled in Matlab). The model of the cell stack is based on Refs. [80][81][82][83] to describe the steady-state and the dynamic stack behavior, respectively. The polarization curve accounts for the activation, ohmic and concentration overpotentials, as reported in Table 5. The open-circuit voltage is evaluated through the Nernst equation, where average partial pressures between inlet and outlet are used; the activation polarization is expressed by using the Butler-Volmer equation, where the charge transfer coefficient is assumed to be equal to 0.5; the ohmic loss is given by the electronic and ionic contribution, R e and R m , respectively; the concentration polarization is evaluated through the hydrogen and oxygen concentration at the electrolyte interface q m . The values of the parameters describing these contributions are based on the references reported in Table 5, and the resulting polarization curve validated against the experimental data reported by Suermann et al. in [56]. A detailed description of the PEME model, considering both steady-state and dynamic behavior, was presented in [42].
The main equations implemented to model the electrochemical devices, as well as the main parameters for process simulations are summarized in Table 5.