Computers and Chemical Engineering Model compendium, data, and optimization benchmarks for sector-coupled energy systems

optimal control of sector-coupled energy systems, there is a lack of available benchmarks. Hence, this article provides a model compendium, exemplary realistic data sets, as well as two case studies (i.e., optimization benchmarks) for an in-dustrial/research campus in an open-source description. The compendium includes stationary, quasi- stationary, and dynamic models for typical components as well as linearization schemes relevant for optimization of design, operation, and control of sector-coupled energy systems.


Introduction
Realistic mathematical models of sector-coupled energy systems are a key for developing tailored numerical optimization methods, which in turn are essential for manifold research efforts towards a successful decarbonization, defossilization, and decentralization of energy supply. Numerical optimization allows to optimally plan, design, operate, and control energy systems while accounting for the inherent volatility of renewables as well as for environmental, economic, and social aspects, see Andiappan (2017) ; Mitsos et al. (2018) .
In process systems engineering, numerical optimization is in many cases the method of choice for control and automation problems ( Engell, 2007;Engell and Harjunkoski, 2012;Kadam and Marquardt, 2007 ). Moreover, its importance for design optimization ( Frangopoulos, 2018 ) and operation of energy systems is steadily increasing. The respective research efforts regarding the dynamic optimization of energy systems comprise a variety of methods and applications, spanning from the development of accurate and fast simulation methods for the control of thermal energy storage ( Barz et al., 2018 ) via the incorporation of real-world weather forecasts ( Constantinescu et al., 2011 ) to nonlinear modelpredictive control and/or real-time optimization of power grids with storage ( Adeodu et al., 2019;Braun et al., 2018;Faulwasser and Engelmann, 2019;Matke et al., 2016 ), and the optimization of HVAC (Heating, Ventilation, and Air Conditioning) systems for buildings ( Bürger et al., 2018;Harb et al., 2015;Perez et al., 2016;Touretzky and Baldea, 2016;Zhang et al., 2014 ).
Design of energy systems is typically cast as Mixed-integer Nonlinear Programs (MINLPs), e.g., in Li and Barton (2015)  Mixed-integer Linear Programs (MILPs), e.g., in Lara et al. (2018) ; Mancarella (2014) ; Zhang et al. (2019) . Besides, dynamic optimization problems comprising operational optimization and optimal control are typically solved by direct methods based on discretization, yielding MINLPs, NLPs (Nonlinear Programs), or MILPs, see Biegler and Grossmann (2004) ; Grossmann and Biegler (2004) . Optimization problems arising in context of sector-coupled energy systems are challenging for a number of reasons: multiple time scales, large number of equations, uncertainties, potentially conflicting multiple objectives, or the hard-to-quantify effect of underlying modeling assumptions. Hence, there are widespread and ongoing research efforts on modeling and numerical optimization for such systems, see, e.g., Barton (2009) ;Faulwasser et al. (2018) ; Majewski et al. (2017) ; Mühlpfordt et al. (2018) ; Roald et al. (2017) .
The development of any numerical method benefits from welldefined, transparent, and realistic benchmark problems. In numerical optimization, benchmark libraries are therefore widely established, including MINLPLib ( Bussieck et al., 2003 ), PrincetonLib ( Vanderbei et al., 2004 ), COCONUT benchmark ( Shcherbina et al., 2002 ), MINTOC benchmark ( Sager, 2012 ), and MIPLIB ( Koch et al., 2017 ). Benchmark problems are also common in Process Systems Engineering such as the Williams-Otto reactor ( Williams and Otto, 1960 ), which, up to this day, is frequently used to compare methods for real-time optimization of process systems ( Srinivasan and Bonvin, 2019 ). Another well-known benchmark problem is the Tennessee Eastman process proposed by Downs and Vogel (1993) . It is still in use for a wide range of research purposes, e.g., for demonstrating the efficiency of a newly developed plantwide control scheme ( Luppi et al., 2018 ). Recently, similar effort s have been made for energy systems. This includes • software frameworks for modeling and optimization of sectorcoupled energy supply systems available online, compare for example Augenstein et al. (2005) , Ringkjøb et al. (2018) , Schütz et al. (2017) , and the Temoa framework with the widely accepted linear and quasi-stationary benchmark model Utopia ( Howells et al., 2011;Hunter et al., 2013 ); • a lately increasing number of open source data-bases providing country-scale data like the Open Energy Modeling Initiative ( Open Energy Modelling Initiative, 2019 ) or specifications of complete power plants like the JRC Open Power Plants Database ( Hidalgo Gonzalez et al., 2019 ); • benchmarks for specific energy systems like electricity grids ( IEEE, 2018;Hörsch et al., 2018;University of Washington, 2018 ) or energy supply systems for supermarkets as used, e.g., in Beykal et al. (2018) ; and • specifically focused model collections, e.g., for ship energy systems ( Sakalis et al., 2019 ) or for economic design of combined cooling, heating, and power (CCHP) systems ( Rech, 2019 ).
However, there are also certain shortcomings and gaps in the benchmarks available in the literature.
• The existing software frameworks usually lack transparency.
Typically, there is no easily accessible documentation of model equations and corresponding example data. • Existing data-sets do not account for the scale of an industrial/research campus or specifications of individual components like boilers. • Specific benchmarks and model collections lack flexibility. If the focus is on a particular sector-coupled energy system, different objectives, or the combination of quasi-stationary and dynamic operation of different components, the adaptation of the proposed setting is usually time consuming or even effectively impossible due to a lack of documentation.
To the best of the authors' knowledge, there is currently no widely accepted benchmark for the optimization of design, oper-ation, and control of energy systems supplying industrial/research campuses. Moreover, models and especially input data for such energy systems are scattered over numerous publications, specification sheets, websites, etc. and are often subject to data protection regulations.
Thus, we compile a model compendium for typical components of an energy supply system coupling cooling, heating, and electricity for industrial/research campuses. The considered setting is inspired by the real-world supply systems of the Campus North of the Karlsruhe Institute of Technology and the Forschungszentrum Jülich. 1 We provide a complete data set for weather, energy prices, cooling, heating, and electricity demands as well as parameter values for the system components. On the demand side, we include building models validated with real-world measurements. The application setting considered is the optimal design and operation of industrial/research campuses. In principle the scope of the models could be extended to include industrial processes and large-scale networks; however, this is outside the scope of the manuscript. Finally, we propose two optimization case studies as open-source benchmark problems: (1) a bi-objective design optimization of a generic energy system and (2) the operational optimization of a dynamic model for a sector-coupled energy system with fixed design.
We address requirements of different application contexts, namely the optimization of design, operation, and control of energy systems. Hence, we include stationary, quasi-stationary, and dynamic model equations of the most frequently considered components. We do not attempt an extensive review on modeling strategies or data ranges. Rather we present one common model formulation for each component. These models are based on already available publications and specification sheets. In fact, the components are modeled based on energy flow rates, which have to satisfy energy balance equations and input-output-relations given by efficiency or COP (coefficient of performance) curves, i.e., they are first-principles models. For the sake of self-containment, we recall a linearization scheme  and an approach for considering a minimum load fraction within the purely continuous and smooth dynamic optimization model of our second benchmark case study based on standard techniques. To foster accessibility of the compendium, we propose a consistent notation.
The remainder of this article is structured as follows: In Section 2 , we present the virtual campus including all components of an energy supply system, which are considered in this article. Based on this generic energy system , we emphasize how the consistent notation and modeling of our corresponding model compendium allows for identifying synergies and structural differences of the various fields of applications. Section 3 contains description and results of the two benchmark case studies. Finally, conclusions are drawn in Section 4 . The Appendix starts with the introduction of the notation in Appendix A . Afterwards, the model compendium including example parameter values forms Appendix B . In Appendix C and Appendix D , the linearization scheme and a smooth extension of the dynamic equations regarding a minimum load fraction, respectively, are given. The git repository available at https://git.es2050.org/heci/energy-benchmark contains input and output data as well as the models and optimization formulations in GAMS ( McCarl and Rosenthal, 2016 ), Modelica ( Mattsson and Elmqvist, 1997 ), and Pyomo http://www.pyomo.org/ . Fig. 1 shows the generic structure of a sector-coupled energy supply system for a campus of variable size considered in this arti-cle. The eventual choice of components, number of respective units as well as coupling points between the gas, electricity, heating, and cooling grid is subject to the specifically chosen application. Both economic and environmental criteria can be chosen for evaluating design and operation of the generic energy system.

Generic energy system for a campus
To cover the specific needs of the optimization of design, operation, and control of energy supply systems, consistent quasistationary and dynamic models as well as linear and nonlinear models are provided. The quasi-stationary models can be obtained by setting all temporal derivatives to zero in the dynamic models. Similarly, the linear models converge to the nonlinear models in the limit of zero distance between two adjacent supporting points. As a result, the gathered component models form a unified framework, which allows to investigate the impact of linearizations and the assumption of quasi-stationary behavior for the optimization of energy supply systems without distortions caused by different modeling premises. Note that we present only models based on energy flow rates. Moreover, we provide dynamic extensions for models with thermal inertia on a scale of minutes. In contrast, models referring to electricity are assumed to be (quasi-)stationary, since the dynamics of electrical processes is usually on a scale of seconds.
The model compendium aims to emphasize synergies and structural differences between optimization applications in energy supply systems. For instance, the incorporation of battery models intrinsically introduces nonsmoothness/nonconvexity in any optimization problem, cf. discussion in Appendix B.5.2 . Moreover, some degrees of freedom of a model depend on the type of component rather than the field of optimization: the input power/heat transfer rate is strictly determined by the output for conversion components, only bounded by the technical possible output for generation components and loosely coupled to the output via storage level and capacity for storage components. However, other degrees of freedom depend on the field of application and increase the complexity of the optimization problem significantly, even if the same model equations are used. E.g., the linearization scheme of the efficiency or COP in Appendix C.1 is sufficient for operational optimization problems with a-priori known capacities, while it needs to be combined with Glover's reformulation ( Glover, 1975 ) in design optimization problems where capacities are degrees of freedom. Besides, the choice between algebraic and dynamic optimization can influence the mathematical properties of the same physical variable. As an example, the roles of input and output heat transfer rate of a boiler are interchangeable for the optimizer within a quasi-stationary, algebraic problem, see Eq. (B.4a) , while their hierarchy is fixed in the dynamic model, see Eq. (B.5) .

Optimization benchmarks
As a proof-of-concept as well as to show the wide range of applications of the provided model compendium, we propose two case studies as optimization benchmarks: (1) a bi-objective design optimization accounting for economic and environmental criteria based on all components of the generic energy system depicted in Fig. 1 , and (2) a dynamic operational optimization of a sectorcoupled energy supply system with fixed, optimal design. Both case studies consider the demand of six office buildings of type "OB", two smaller office buildings of type "OBM", and two experimental facilities "EF", each with one thermal zone, see Appendix B.2.2 for details. More specifically, Appendix B.2.1 describes the parameter identification of the thermal models for the office building "OB 1" as well as the experimental facility "EF 1" based on measured temperature data. Figs. 2 and 3 show an excerpt of the simulated temperature and heating/cooling input Q dem z compared to the real measurement data for buildings "OB 1" and "EF 1". The Coefficients of Variation of Root Mean Square Error CV(RMSE) (acc. to Coakley et al. (2014) ) of the indoor temperature are 1.2% for "EF1" (for hourly data between 13 Jan and 11 Sep 2018) and 0.3% for "OB1" (for hourly data between 1 Jun and 31 Dec 2018).
Example values for general parameters and inputs as well as all model equations can be found in Appendix B . Parameter values which differ for the two case studies are given in the following. This includes, e.g., bounds on capacities which are required for the design optimization and fixed to the nominal value for the optimization of the operation.
We scale the values of system variables within the optimization models to a range of approximately 0 to 1 to avoid numerical problems.

Design optimization
Design optimization has to cope with large-scale optimization problems, in particular, due to the incorporation of combinatorial decisions and operational optimization ( Frangopoulos et al., 2002 ). Goderbauer et al. (2019) have even shown that the design problem of (distributed) energy supply systems is NP-hard. In this case study, we regard both minimum costs and minimum global warming impact, which further increases the complexity of the design optimization problem. Thus, the models are linearized as described in Appendix C , allowing for continuous sizing of all components. Additionally, the time-varying input parameters are aggregated using the k-medoids method proposed by Bahl et al. (2018b) . Thereby, we employ 4 typical periods and 4 segments per typical period with additional peak values for demands. The resulting sorted aggregated time-varying demands and prices are shown in Fig. 4 and are available at https://git.es2050.org/heci/ energy-benchmark in directory "3_1_Design_Optimization" in file "AggregatedTimeSeries.csv". Despite the clear deviations between the aggregated and the original full time-series, aggregating timeseries have been shown to lead to near-optimal solutions in studies ( Bahl et al., 2018a;. All components of the generic energy system depicted in Fig. 1 are considered for the supply of a campus comprising ten buildings. Only wind turbines are excluded as the installation is often prohibited due to construction limits as, e.g., minimal distances to neighboring residential buildings. The values of Voll (2013) are used for parameters regarding boilers BOI , combined heat and power engines CHP , absorption chillers AC , and compression chillers CC , while the parameter values of  are taken for the photovoltaic units PV , heat pumps HP , batteries BAT , and thermal energy storage units T ES , see Tables 1 and 2 . Note that the corresponding references had demonstrated the linearization to be an adequate representation of the nonlinear models. Example values for linearized part-load behavior and linearized investment costs are given in Tables C.1 and C.2 , respectively. Please note that cyclic conditions apply for the operation of the thermal storage units as well as the battery in each typical period. Further, charge and discharge of all storage units are suppressed in the time steps representing the additional peak demand values. We assume a constant global warming impact for the electricity mix of the grid gwi el = 561 g CO 2 −eq . / kWh . As an alternative, a time dependent global warming impact based on  is provided for the considered time-series as well, see "Aggregat-edTimeSeries.csv" at https://git.es2050.org/heci/energy-benchmark in "3_1_Design_Optimization". When electricity is fed into the grid, we assume a credit for the global warming impact GWI (B.1.2) following the idea of the avoided burden ( Baumann and Tillman, 2014 ). Moreover, we employ gwi fuel = 244 g CO 2 −eq . / kWh for the specific global warming impact of purchased gas as well as time horizon τ h = 4 a and interest rate γ 5 = 8 % for the calculation of the present value factor PVF ( Majewski et al., 2017 ). The interest rate in application cases strongly varies; other authors employ 5 % for example ( Schütz et al., 2017 ).
To solve the design optimization problem, we apply the automated super-structure-generation approach from Voll et al. (2013) . This approach successively performs superstructure optimization until the objective function value does not further improve. In each optimization, the superstructure of the energy supply system is enlarged by one unit for each component type. Although we consider three sizes of combined heat and power engines CHP , see Table 1 , we only allow one additional CHP unit rather than one for each size when increasing the superstructure by one unit. To decrease the computational effort, we aggregate the roof area of the office buildings and of the experimental facilities, respectively, and we allow at most one photovoltaic unit and one storage component for each building type. Note that the values for the linearized investment costs are not changed despite the aggregation, since there is no significant economy of scale for photovoltaic components.
The design optimization approach is implemented in GAMS 24.7.3 ( McCarl and Rosenthal, 2016 ). We choose GAMS as it is one of the standard modeling environments. To solve the problem, CPLEX 12.6.3.0 ( IBM Corporation, 2015 ) is used employing an optimality gap of 0.5 %. A GAMS file containing the problem statement after automated super-structure generation and using the single objective of total annualized costs, namely "Model.lst", can be found at https://git.es2050.org/heci/energy-benchmark in "3_1_Design_Optimization". Therein, we also provide a corresponding pyomo file generated automatically via the GAMS convert function, along with the variable mapping. Pyomo has the advantage that it is open source.
We perform a bi-criteria optimization, minimizing the total annualized costs TAC and the global warming impact GWI employing the augmented ε-constraint method ( Mavrotas, 2009 ). The re-sulting trade-off in the Pareto front as well as the corresponding Pareto-efficient designs are shown in Fig. 5 .
The global warming impact GWI decreases for designs with a tri-generation system in place of separately operating boilers BOI and compression chillers CC with additional purchase of electricity from the grid. Low-GWI energy systems employ a higher number of CHP engines in combination with the installation of absorption chillers AC , which allows for simultaneous heating and cooling supply while providing electricity for on-site demands or the grid. As a result, the system is even able to reach a negative global warming impact GWI , as the avoided burden by feeding in electricity ( Appendix B.1.2 ) is higher than the global warming impact induced by the consumption of fuel on-site.
The increasing investment in a higher number of smaller conversion units and larger storage units enables a more ecological operation by utilizing highly-efficient operational points and load shifting, compare Fig. 5 (b). Moreover, larger photovoltaic units are installed for a more climate-friendly energy supply. The maximum PV area is reached at the third point on the Pareto front with decreasing global warming impact GWI .
Herein, we do not consider the selection of a final design by the decision maker. This selection can be done with a wide range of decision supporting tools, see, e.g., Jing et al. (2019) . For the synthesis of distributed energy supply systems and other two-stage optimization problems, for instance, the flex-hand approach automatically selects a highly flexible design in operation such that the final design performs well regarding all considered criteria ( Hollermann et al., 2019 ).

Operational optimization
The second case study pertains to the offline optimal control of nominal operation based on realistic simulation data and parameter values as well as the dynamic models given in Appendix B with the extension discussed in Appendix D . The energy supply system considered is given by the cost-optimal solution of the design optimization problem discussed in Section 3.1 . We particularly focus on the parallels between thermal energy storage (TES) and the thermal inertia in buildings.
The feasibility of the design for an energy system with either explicit or implicit storage is guaranteed by excluding TES units in the design optimization. Instead, the size of the TES unit is adapted to heat transfer capacity Q max = 100 kW of Building "OB  The controls are the heat transfer input of the boiler and the CHP unit, the power provided by the PV components clustered in one unit for all (modified) office buildings and one for both ex-perimental facilities, as well as the purchased electricity at any considered time point. The initial values of the load fractions of boiler and CHP are optimized. The TES is set to be half-full at the beginning to allow for both charging and discharging. The heat transfer rate exchanged with the TES is the difference between the given demand and the heat transfer rate provided by the boiler and the CHP unit. We minimize the total costs subject to the model constraints reported in Appendix B . Note that violations of path constraints may occur between discretization points, compare Fu et al. (2015) . Besides, we do not impose periodic boundary conditions to not further restrict the optimizer. Since weather is not week-periodic, periodic operation is not expected to be optimal. Other researchers, e.g., ( Ghobeity and Mitsos, 2012 ) have imposed periodic boundary conditions to avoid discharging the storage. Our model compendium enables variations of the benchmark to account for such boundary conditions. Note however that periodic boundary conditions and oscillating systems bring substantial challenges ( Wilkins et al., 2009 ).
For the dynamic optimization of the operation, the dynamic optimizer DyOS ( Caspari et al., 2019;DyOS, 2019 ) calling local nonlinear optimizer SNOPT ( Gill et al., 2005 ) and integrator IDAS ( Serban et al., 2018 ) is employed. We use feasibility tolerance of 0.01 and optimality tolerance of 0.001. We write the model in Modelica. The motivation is that it is open source and supported by a variety of commercial and open-source simulation and optimization tools, including our in-house solver DyOS. In nonconvex dynamic optimization problems, a good initial guess is typically required for convergence of the optimizer. While in principle deterministic global methods exist since more than a decade ( Singer and Barton, 2006 ), these are far from being applicable to such systems. An alternative are heuristic local methods such as multistart, which has been applied to energy systems, e.g., ( Ghobeity and Mitsos, 2012 ). However, a challenge is that many runs fail to converge. Herein, we find the initial point based on ad-hoc adaptations of intermediate optimization results. Flat Modelica files containing the dynamic optimization problem as well as the dynamic simulation model for building "OB 1" are given by "OptModel.mo" and "OB1SimModel.mo", respectively, at https://git.es2050.org/heci/energy-benchmark in directory "3_2_Optimal_Control".
The CHP over-fulfills the electricity demand when the gas price is low and the gas demand is high compared to the electricity demand, cf. Fig. 7 (b). Thus, no purchase of electricity is required and only the volatility of the sale price of electricity can affect the optimal solution. In this way, the TES enables a higher CHP load in periods with higher sale prices for electricity. For instance Fig. 7 (a) shows the CHP extended full-load period within interval [18h, 24h] and the load fraction peak within interval [150h, 168h]. In contrast to that, the boiler is only used during the high demand periods of the workdays, i.e., from 0h to approximately 120h. To satisfy the remaining heat demand, the boiler is run at about 30 % to 40 % load, see Fig. 7 c. Note that the operation with this low load fraction is already highly efficient, cf.  Finally, we transfer the capability of the TES as explicit storage to Building "OB 1" as implicit storage. In particular, the optimized heat transfer rate interchanged with the TES Q i , i ∈ T ES is added to heating demand Q heat,dem given in B.2.2 which has been used for the dynamic optimization. The internal temperature within Building "OB 1" is simulated based on heat demand Q heat,dem , i.e., for an energy system with an explicit storage only, and the combined demand Q heat , dem + Q i , i ∈ T ES , i.e., for an analogous energy system with an implicit storage only. Fig. 8 shows that the transfer of the heat transfer rate Q i , i ∈ T ES from the TES to Building "OB 1" leads to visible but for human hardly sensible oscillations in the range of up to 0.2K.

Conclusion
This article presents a model compendium for common components of energy supply systems present in industrial or research campus areas. Moreover, the included validated building models rely on real-world data from the Campus North of the Karlsruhe Institute of Technology and from the Forschungszentrum Jülich ( HECI, 2019 ). We provide one model formulation on the scale of energy flow rates for each component considered. The model compendium is structured in terms of notation and modeling principles such that it can be extended by additional components, e.g., solar-thermal collectors and power-to-X technologies, or by including high-fidelity models, e.g., for gas grids and thermal grids.
The compendium addresses requirements of different fields of applications, namely the optimization of design, operation, and control of energy supply systems. Hence, it includes both quasistationary and dynamic models as well as linearization schemes. The notation of all given models is unified for more transparency concerning synergies and structural differences of the different fields of applications. This way, we aim to support the transfer of models and methods between the different fields of applications. Moreover, the unified modeling framework allows investigating the influence of different formulations on computation times and on the accuracy of solutions. For instance, two corresponding optimization problems can be obtained for the same individual energy system by replacing the quasi-stationary model of one component by the respective dynamic model.
Additionally, we propose two optimization benchmarks exploiting the wide range of presented model formulations. In the first case study, a bi-criteria design optimization regarding total annual costs and global warming impact is performed for the generic energy system based on linearized quasi-stationary models. The results of the design optimization hint at the benefit of photovoltaic components, storage systems as well as the synergy of trigeneration for an ecological energy supply. In the second case study, the operational optimization of an energy supply system based on nonlinear dynamic models emphasizes the possibility to exploit varying electricity prices with the help of a combined heat and power engine and a thermal energy storage. For the sake of illustration, we also touch upon the role of thermal inertia of buildings via subsequent simulations. The case studies constitute substantial numerical challenges, e.g., for testing global solution methods for the operational optimization. Both can be easily adapted, e.g., to allow for different boundary conditions for the operation.
Notably, the benchmarks come with a complete set of readyto-use input data and the respective model files, namely a GAMS ( McCarl and Rosenthal, 2016 ) listing file for the design optimization and a Modelica ( Mattsson and Elmqvist, 1997 ) file for the dynamic optimization of operation, available at https://git.es2050.org/ heci/energy-benchmark . We also provide equivalent Pyomo files. The data sets may be extended by real-world measurements of demands and their corresponding price and weather data to account for model-plant-mismatch or real-world uncertainties.
The novelty of our approach is the definition of suitable benchmarks, writing consistent models for important unit operations allowing for various use cases, and combining these models with useful data. We utilize established solution methods and the models are not fundamentally different from existing state-of-the-art models. We envision the modular compendium, the nominal data set, and the benchmarks to enable transparent comparisons of optimization methods for sector-coupled energy systems.

Author Contribution
The author contribution possibility was offered after our submission.

Declaration of Competing Interest
We have no conflict of interest.

Acknowledgment
Dedication: This article is dedicated to Lotte who chose a dreamlike world.
Funding: This work was supported by the Helmholtz Association under the Joint Initiative "Energy System 2050 -A Contribution of the Research Field Energy". Susanne Sass is grateful for her association to the International Research Training Group (DFG) IRTG-2379 Modern Inverse Problems.
We are thankful to John Siirola for his helpful suggestions on Pyomo and to Bethany Nicholson for her help in converting Modelica to pyomo.dae.
Author Contributions: SS led and coordinated the writing to which all authors contributed. SS, TF, DEH, DS, and TS equally contributed to the model compendium. AB, LG, VH, DM, and AM guided the research of SS, DEH, DS, and TS. DS performed the parameter identification for the models of buildings and the corresponding simulation of demand data with input and guidance from LG. DEH performed case study "Design optimization" with input and guidance from AB. DYS prepared the publishable material. SS performed case study "Operational optimization" with input and guidance from AM. CK prepared the publishable material.

Appendix A. Notation
The variables, see Table A.1 , are specified with the help of superscripts, see Table A.3 , and assigned to certain sets via subscripts, see Table A.2 . Apart from this, we omit explicit time-dependency in equations for the sake of simplicity unless it may create confusions. Thus, variables which depend explicitly on time are written in standard font, e.g., load fraction λ, while all other variables are represented by bold symbols, e.g., minimum load fraction λ min or efficiency η(λ) . Whenever time-dependent variables occur in equations, these equations have to be satisfied at any considered time point t ≥ 0.

Appendix B. Model compendium
In this section, economic and environmental evaluation criteria for design, operation, and control of a sector-coupled energy supply system as well as (non-)linear quasi-stationary and dynamic models are given for common components. This includes models for office buildings and experimental facilities; models for the conversion components boiler ( BOI ), combined heat and power engine ( CHP ), absorption chiller ( AC ), turbo-driven compression chiller ( CC ), and heat pump ( HP ); models for photovoltaic units ( PV ) and wind turbines ( W T ) clustered as generation components; models for thermal energy storage ( T ES ) as well as batteries ( BAT ); and, finally, models for the thermal and the electricity grid coupling these components.
The linearization scheme for the nonlinear efficiency curves and investment costs is given in Appendix C . Finally, an extension of the dynamic model of boiler and CHP for the inclusion of a minimum load fraction without the introduction of binaries or nonsmoothness is given in Appendix D .

B1. Evaluation criteria
In the following Sections B.1.1 and B.1.2 we introduce the total annualized costs TAC and the global warming impact GWI as evaluation criteria.

B1.1. Total annualized costs
Total annualized costs (TAC) are an economic criterion for the evaluation of an energy supply system with respect to operational, investment, and maintenance costs. The TAC can be calculated by with duration of time step 8760 h · τ dur per year τ dur , prices c fuel , c el,buy , c el,sell , purchased energy rate with respect to natural gas Q fuel,in , purchased and sold electricity P buy and P sell , respectively, present value factor PVF , factor for maintenance costs per year γ 4, i , and capital expenditure CAPEX i .
The present value factor can be calculated by ( Broverman, 2010 ) PVF = ( γ 5 + 1) τ h − 1 ( γ 5 + 1) τ h · γ 5 with interest rate γ 5 , e. g., γ 5 = 8 % , and time horizon τ h , e. g., τ h = 4 a . In this study, we obtain the CAPEX by the power law of capacity Smith (2005)   Instead of the total annualized costs TAC , any other economic objective function could be chosen. However, when regarding economic and ecologic criteria, the total annualized costs TAC lead to Pareto fronts with low environmental impacts compared to other economic objective functions ( Pintari č and Kravanja, 2015 ).

B1.2. Global warming impact
An environmental evaluation criterion of energy supply systems is given by its global warming impact where gwi fuel and gwi el are the specific global warming impacts of the energy sources, for example values see Section 3.1 or Federal Environment Office (2018) . Note that the specific global warming impact of purchased electricity is varying remarkably over time. We follow the idea of the avoided burden ( Baumann and Tillman, 2014 ) and assume a credit for the global warming impact GWI when electricity is fed into the grid. As the operation usually affects the global warming impact significantly higher than the manufacturing of the components ( Guillén-Gosálbez, 2011 ), we only consider the contribution of the operation.

B2. Buildings
The models of buildings origin from the energy balance equation based on the following assumptions: (A1) Only internal energy is considered. In particular, the kinetic energy of the system is neglected and potential energy cancels out. (A2) The change of internal energy equals the heat flow, i.e., no additional work is applied to the system. (A3) The mass balance is fulfilled.
As sketched in Fig. B.1 , the physical effects which most strongly influence the temperature within the buildings are identified as (P1) Heat transport mechanisms with external air of temperature T amb via air exchange based on air change rate γ amb , k 7 and heat capacity C air, z as well as heat transfer through walls, windows, roof, and floor based on heat transfer coefficient γ amb , U 7 , (P2) Heat input by solar irradiance Q irr s,z on the building surface s with the solar energy absorption coefficient γ irr 8 according to Harb et al. (2016) , and (P3) Installed heating/cooling system with heating/cooling input Q dem z with heating factor γ th 9 . (P4) Heat capacity C z takes into account internal walls, air, external walls, roof, and basement floor. However, it is only a fraction of the sum of all heat capacities C tot of the components of the building. This is presumably the case, since the outer shell is stronger coupled to the ambient temperature than to the inside and therefore does not significantly contribute to the indoor climate. Due to missing data, the Eq. (B.3) provides a dynamic equation for differential state T z in dependence of control input Q dem z as well as time-varying parameters T amb and Q irr s,z , s ∈ S for any thermal zone z ∈ Z. Example values are given in Section B.2.1 .
The model is based on Park et al. (2011) . It is chosen to be as simple as reasonable for easy integration in different use cases and has been extended by considering solar irradiance. The parameter values given serve as first orientation for researchers without access to building models. It is generally advisable to estimate the parameters on data as the parameters vary depending on the characteristics of the buildings. The model can be adopted to different building sizes or orientations by changing the surface areas and the dependent solar irradiance. Moreover, in the case studies, we consider simple cuboid shape. If the shape of the building is different, it would be sensible to estimate the parameters to temperature measurements and/or explicitly take self-shading effects into consideration.
Note that parameter identification via linear regression of averaged measurement data of the indoor temperature shows that the consideration of additional factors like human body heat, electrical devices, wind velocity, coupling between different zones, and wall temperatures leads to worse identification results due to linear dependencies displayed by high condition numbers.

B2.1. Parameter identification
The parameters ( γ amb 7 ) z , γ irr 8 , and C z of the gray-box model are fitted to measured temperature data T z with γ th 9 ≤ 1 and the input time series T amb , Q dem z , and Q irr s,z . In fact, dimensions and common information of the two original buildings, the office building "OB 1" and the experimental facility "EF 1", are retrieved from internal reports as shown in Fig. B.2 and Table B.2 .
Moreover, ambient temperature T amb is given by weather data of year 2018 from Deutscher Wetterdienst for Stuttgart ( DWD Climate Data Center (CDC), 2018 ). Solar irradiance I s is calculated using the horizontal global irradiance and the horizontal diffuse irradiance of the weather data as well as considering the orientation  of each surface s ∈ S by angle φ s according to Kreider et al. (2010) .
The two original buildings have no solar panels installed, therefore possible shading effects are neglected. Finally, solar heat input Q irr s,z is the product of solar irradiance I s and the respective area A s,z for any surface s ∈ S and thermal zone z ∈ Z. The values resulting from the parameter identification process with one thermal zone are given in Table B.3 . Note that the parameters for buildings "OB 1" and "EF 1" are identified based on real measurement data. For the parameters of "OB 2" to "OB 6" as well as "OBM 1" and "OBM 2", intervals are predefined according to empirical considerations. Their parameters approximately follow a uniform distribution within those intervals. The parameters of "EF 2" are slight variations of those of "EF 1". The following Section B.2.2 depicts the simulation of the heating/cooling demand.

B2.2. Simulation of building demands
For simulating realistic building demands 3 , the common approach of a standard PI controller representing a thermostat is used for the control of the buildings ( Peeters et al., 2008 ). The deviation from desired temperature T 0 is the input. The feedback loop for any thermal zone z ∈ Z is defined by 3 These demands are often referred to as building loads in the community for modeling, simulation, and optimization of buildings.    The aggregated electrical demand profile for the ten buildings is the sum of 1) three generated demand profiles following the G0 demand profile of BDEW (2019) with Saturdays treated like Sundays as well as different offsets and gains for each demand profile for three OB/OBM buildings; 2) measured data of "OB 1" for the years 2014 to 2018 which are shifted to start at the same day of week in order to retrieve five more OB/OBM demand profiles; and 3) measured data of "EF 1" for the years 2017 and 2018 to get two demand profiles for the EF buildings.

B3. Conversion components
The quasi-stationary model for general conversion components C = BOI ∪ CHP ∪ HP ∪ AC ∪ CC is given as in Voll et al. (2013) by (B.5) see Sass and Mitsos (2019) for more details. Note that this dynamic model is based on a simplified energy balance that considers heat transfer rates rather than temperatures. Furthermore, all heat losses are assumed to be proportional to the input heat transfer rate.
The efficiency or COP curves η of the respective components are given in Table B.5 . Note that the given curves are used for both the quasi-stationary and the dynamic models. Apart from this, all efficiency and COP curves are assumed to be temperature independent, which is reasonable for boilers and CHPs but not necessarily for chillers and heat pumps ( Augenstein et al., 2005 ).

B4. Generation components
In the context of this article, the output of generation components is limited by their capacity and the availability of renewable energy resources, namely solar irradiation and wind. However, the maximum available power may not be exploited, e.g., if this would impair grid stability or exceed storage capabilities.

B4.1. Photovoltaic units
The electrical power P i provided by a photovoltaic (PV) unit i ∈ PV is limited by the solar irradiance I , the total area A i of the unit and its efficiency η i via Thereby, I accounts for direct, diffuse, and reflected solar irradiance onto the tilted PV unit's area.
Furthermore, the PV unit cannot exceed its nominal capacity where the nominal capacity depends on the area of the unit as well ( Ren et al., 2009 ). The average and maximum efficiency of current German PV technologies are 17% and larger than 20%, respectively, while the average performance ratio ranges from 80 to 90% ( Wirth and Schneider, 2019 ). Thus, we choose efficiency η i = 0 . 19 and nominal capacity P nom i = 0 . 171 kWm −2 · A i as an example.
For the case studies, the PV components are virtually installed on the roofs of the buildings described in B.2 with the parameters given in Table B.2 . To avoid self-shading effects of the panels, it is assumed that the rows have a minimum distance of three times the (projected) height of the units. Thus, for an inclination of 30 • the maximum available PV surface area is 42% of the total roof area, provided that the complete width of the roof can be used to install the panels. Applying that rule to the two experimental facility buildings, we have A max EF ≈ 205 m 2 . The PV area for the office buildings is up to 85% of the total roof surface, i.e. A max OB ∪ OBM ≈ 5261 m 2 , since the inclination is only 10 • . The irradiance I is calculated as described in Section B.2.1 . The data is included in WeatherAndDemandTimeSeries.csv available at https://git.es2050.org/heci/energy-benchmark in directory "App_Weather_And_Demand".

B4.2. Wind turbine
The maximum power output of a wind turbine is determined by the wind velocity, which corresponds to the part-load behavior of the wind turbine, and its nominal capacity. We introduce efficiency η el W T (λ i ) for the mapping of wind velocity to the power output for each wind turbine i ∈ W T and, thus, obtain As an example, the efficiency given by ENERCON (2015) may be used. Thereby, load fraction λ i is the wind velocity at each time step, normalized by the rated wind velocity, e.g. 15 m/s.

B5.1. TES
In this article, we only take into account a lumped model of a hot/cold water storage tank for the thermal energy storage (TES). More sophisticated multi-layer tank models are discussed in, e.g., Steen et al. (2015) ; Schütz et al. (2015) .
The energy balance based on heat transfer rates of such a simple TES model yields with added and withdrawn heat transfer rates Q in i and Q out i , respectively, efficiencies η in i and η out i as well as self-discharge in dependence of the currently stored energy E i with time constant With reformulation Eq. (B.8) the number of degrees of freedom is reduced, since only the total heat transfer rate flowing through the thermal energy storage is considered. Note that the dynamics in Eq. (B.7) are commonly discretized using the implicit ( Schütz et al., 2017 ) or explicit Euler scheme ( Bahl et al., 2018b ). However, in benchmark case study "Operational optimization" we stick to formulation Eq. (B.8) , since a more sophisticated integration scheme is incorporated in the dynamic optimization framework used. Aside from that, binary variables can be introduced to prevent simultaneous charging and discharging, cf. battery model in Section B.5.2 . The storage tank's capacity E nom i serves as an upper bound for the TES and presents a design variable determining its capacity In the design optimization, capacity E nom i is limited by (B.10) Moreover, the heat transfer rates for charging and discharging the TES are limited by (B.12) with rates 1/ τ in and 1/ τ out limiting the charging and discharging process, respectively. Appropriate values are given by 1 / τ in = 1 / τ out = 1 h −1 , cf. ( Bahl et al., 2018b ).

B5.2. Battery
A generic model of an electrical battery is given by replacing the heat transfer rates in Eq. (B.7) by the electrical power applied to the battery. However, self-discharge is often negligible for (Li-Ion) batteries ( Zimmerman, 2004 ). This yields for any battery i ∈ BAT Constraints Eqs. (B.9) to (B.12) apply analogously to the battery model if heat transfer rate Q is swapped with power P where applicable. In contrast to a TES, the charging and discharging process of batteries is typically not limited by rate constraints.  report time constants τ in = τ out = 4 . 2 · 10 −5 h , which lead to large upper bounds in Eqs. (B.11) and (B.12) . Note that the constraints implied by underlying power electronics are usually considerably tighter.
Similar to the case of TES models, cf. Section B.5.1 , it is quite common to consider the discrete-time counter part of Eq. (B.13) in scheduling of power systems. To this end, the ODE can be discretized by the forward Euler method considering averaged values of P in ( t ) and P out ( t ) and a constant step width of e.g. 15min. We remark that the given model does not account for specifics of all existing battery technologies. For example, detailed models for RedOx-Flow ( Blanc and Rufer, 2008 ) and other battery types are beyond the scope of this work.
In contrast to TES models, a battery cannot be charged and discharged at the same time. Put differently, the battery cannot actively dissipate energy. Hence, the constraint is added. As this constraint leads to feasible sets with nondifferentiable boundaries, it has been suggested to either neglect it ( Braun et al., 2018 ), or to relax it as follows The integer variable z i ∈ {−1 , 1 } allows to discriminate charging and discharging. Together with aggregating input and output, the following mixed-integer formulation dE i dt

B6. Grid models
We consider thermal, electricity, and gas grids as components and use simple models. In particular, gas is only a potential energy resource and the gas grid is approximated as a point source with a known gas price.

B6.1. Thermal grid
We do not consider an external thermal grid. Thus, cooling and heating supply have to match the aggregation of building demands and storage capacities.
According to Mehleri et al. (2012) , an energy balance is formulated for each node j of the thermal grid, comprising generation, consumption, storage, and interaction with neighboring nodes l . For a given energy supply system, this results in for each node j of the heating grid and for each node j of the cooling grid, where parameter γ is a loss factor. Mehleri et al. (2012) and Obara (2007) further approximate the loss factor γ proportional to the distance to the neighboring node l .

B6.2. Electricity grid
A simplified model of a balanced electrical AC (alternating current) grid can be given by a lumped-parameter system at steadystate, which can be described by the triple (N , G, Y ) , where N = { 1 , . . . , N} is the set of buses (nodes), G ⊆ N is the non-empty set of generators, and Y = G + j B ∈ C N×N is the bus admittance matrix ( Grainger and Stevenson, 1994 ). The off-diagonal entries of Y can be written as y li = g li + j b li , whereby g li is the conductance for the line li , respectively, b li is the line susceptance. The diagonal entries of Y are y ll = y l + l = m y li , where y l accounts for linear load connected to bus l .
For the sake of simplicity, we assume that there is only one generator per bus (i.e. G ⊆ N ). Thus, at each bus i ∈ N we have The parameter P dem i models the demand of electrical power at node i , it also captures uncontrollable renewable generation, e.g., the maximum power output of PV components. Batteries are considered to be generators.
To reduce nonconvexities, lossless lines, small phase differences, and constant voltage magnitudes are commonly assumed. With these assumptions, the overall active power balance for the grid reads Note that power balance Eq. (B.15) is used in the given benchmark case studies.
A simple expression for the phase angles θ i at each bus is given (B.16) where B is the imaginary part of the bus admittance matrix Y, P is the vector of electrical powers, and θ is the vector of phase angles.
The above Eqs. (B.16) are the so-called DC (direct current) power flow equations ( Grainger and Stevenson, 1994 ).

C1. Linearization of efficiency/COP curves
In this article, we assume that efficiency or COP η i is a given, possibly nonlinear function in dependence of load fraction λ i := Table B.5 . According to Voll et al. (2013) , it is favorable to linearize the functional dependency of input heat transfer rate Q in i on output heat transfer rate Q out In this way, at most one interval can be active j∈J λ b λ i, j ≤ 1 ∀ i ∈ U and the resulting load fraction λ i is obtained by Based on supporting points λ in i, j , λ out i, j , the slope parameter β λ i, j of the line segments j is given by This yields the piecewise linear formulation (C.2) Note that load fraction l i,j equals 0 if the binaries b λ i, j of all intervals ∀ j ∈ J λ are 0. Example values for the linearization of the functions given in Table B.5 is given in (C.4)  ( Glover, 1975 ). Thus, we substitute the bilinear product in Eqs. (C.3) and (C.4) with a new time-dependent variable d i,j and obtain On the one hand, Eq. (C.5) guarantees that the new variable d i,j is 0 if the corresponding binary b λ i, j equals 0. On the other hand, variable d i,j equals capacity Q nom i if binary b λ i, j equals 1 due to Eq. (C.6) .

C2. Linearization of investment costs
Analogously to C.1 , we can replace the nonlinear investment costs , j ∈ J Q for any unit i ∈ U. Example values of the supporting points CAPEX lb i, j , Q lb i, j for the considered components are given in Table C.2 . Note that the nonlinear function CAPEX i only depends on the independent variable Q nom i apart from constant parameter values and, thus, is either constant (operational optimization) or univariate (design optimization). In contrast to C.1 , the application of Glover's linearization is therefore not necessary for the linearization of the investment costs CAPEX i .

Appendix D. Smoothing of efficiency regarding minimum load fraction
In this section, we include the minimum load fraction into the dynamic model equations of boiler and CHP, see Section B.3 , to show the expandability of the presented model equations. More precisely, the hyperbolic tangent is used as a smooth switching function between turned off mode and an operation with a load fraction larger than the minimum load fraction. Thus, no binary controls are introduced and, in particular, we can avoid adding a combinatorial complexity to a dynamic model. , i ∈ BOI ∪ CHP as given in Table B.5 without a binary control turning the boiler or CHP explicitly on. Therefore, the zero-operation below the minimum load fraction is replaced by a linear operation with sufficient slope, see