Demand Response Analysis Framework (DRAF): An Open-Source Multi-Objective Decision Support Tool for Decarbonizing Local Multi-Energy Systems

: A major barrier to investments in clean and future-proof energy technologies of local multi-energy systems (L-MESs) is the lack of knowledge about their impacts on proﬁtability and carbon footprints due to their complex techno-economic interactions. To reduce this problem, decision support tools should integrate various forms of decarbonization measures. This paper proposes the Demand Response Analysis Framework (DRAF), a new open-source Python decision support tool that integrally optimizes the design and operation of energy technologies considering demand-side ﬂexibility, electriﬁcation, and renewable energy sources. It quantiﬁes decarbonization and cost reduction potential using multi-objective mixed-integer linear programming and provides decision-makers of L-MESs with optimal scenarios regarding costs, emissions, or Pareto efﬁciency. DRAF supports all steps of the energy system optimization process from time series analysis to interactive plotting and data export. It comes with several component templates that allow a quick start without limiting the modeling possibilities thanks to a generic model generator. Other key features are the access and preparation of time series, such as dynamic carbon emission factors or wholesale electricity prices; and the generation, handling, and parallel computing of scenarios. We demonstrate DRAF’s capabilities through three case studies on (1) the DR of industrial production processes, (2) the design of battery and photovoltaic systems, and (3) the design optimization and DR of distributed thermal energy resources. are considered within model generation through topology sort. This ensures that all submodels contributing to a collector are built before that collector is executed.


Introduction
Following the Paris Agreement 2015 [1] and the Glasgow Climate Pact 2021 [2] there is an urgent need for decarbonization worldwide. With the European Green Deal approved in 2020, the European Union (EU) aims to achieve carbon neutrality by 2050 [3]. Energy system modeling is essential to driving the rapid adoption of smart energy services and clean innovative technologies needed to decarbonize energy systems [4]. Powerful energy system modeling frameworks have been presented as proprietary and open-source software. Most open-source frameworks, though, were developed with large-scale energy systems in mind, and only a few focus on decision support for local multi-energy systems (L-MESs). However, significant adoption of clean technologies and demand response (DR) programs are required in individual L-MESs, as the success of the energy transition depends on decentralization. According to the International Energy Agency, 70% of clean energy investments over the next decade will have to be made by private developers, consumers, and financiers [5]. Industrial and commercial electricity consumers can contribute to the expansion of renewable energy sources (RES), e.g., through on-site photovoltaic systems. At the same time, they can help to integrate the RES through the smart use of flexible loads and energy storage through DR, which additionally offers cost-saving potential for companies. This demand-side flexibility is needed, as the flexibility demand of the electricity system is expected to quadruple by 2050, even though the availability of conventional flexibility sources will decrease due to the decommissioning of fossil power plants [6]. Therefore, in this paper, we focus on how cost and emission reductions that arise from applying DR to L-MESs, especially in the industrial and commercial sectors, can be quantified using the Demand Response Analysis Framework (DRAF).
Decarbonization of L-MESs are often considered by minimizing greenhouse gas emissions and costs within multi-objective modeling [7]. E.g., in [8], an efficient energy management model based on a multi-objective optimal power flow problem is proposed that considers flexibility of storage units of industrial networks. In [9], the contributions of energy storage systems, production buffer stocks, and smart transformers to a net zero energy factory were analyzed. Andiappan [10] reviewed mathematical optimization approaches for energy system synthesis and identified the concept of eco-industrial parks to be a future research direction.

Demand Response
Due to temporal variability, uncertainty, and location constraints, the integration of high shares of RES is demanding and requires high operational flexibility of the power system, which could lead to integration costs of RES [11]. In DR, final consumers provide demand-side flexibility to the electricity system by voluntarily changing their load profiles in reaction to price signals (price-based DR) or specific requests (incentive-based DR) [12]. It is widely acknowledged that DR is a key element of the transformation to a carbon-free energy system enabling cost-efficient integration of fluctuating RES [13]. Gils [14] identified a significant DR potential in Europe; the minimum aggregated hourly averages of load reduction and increase were 61 and 68 GW, respectively. For a general overview of demand response, the reader is referred to [15].

DR in the Industrial and Commercial Sector
In recent years, several research projects investigating the DR potentials in the industrial and commercial sectors have been funded at the European and national level. Examples on European level are: demand response in industrial production (DRIP) [16], demand response integration technologies (DRIvE) [17], using the flexibility potential in energy intensive industries to facilitate further grid integration of variable renewable energy sources (IndustRE) [18]. Examples on national level are: SynErgie [19], Pilot Project DSM Bavaria [20], refrigerated warehouses store energy for smart energy grid (FlexLast) [21]. Commercial and industrial enterprises have large DR potentials [22], which can be divided into cross-sectional technologies [23] and energy-intensive production processes [24]. However, in order to exploit the full environmental and economic DR potential in the industrial and commercial sectors, aspects such as dynamic pricing [25], sector coupling, onsite generation, and onsite energy storage need to be considered in an integral analysis. Due to this complexity, there is still a lack of knowledge about existing flexibility, which is a major barrier to participation in DR programs [26]. To address this complexity, mathematical optimization, which has already been long proven in the analysis of large international power systems, can be applied to distributed energy systems. However, data preparation, model formulation, scenario definition, and result presentation require relevant experience and expertise to which the decision-makers in the companies often do not have access.
The aim of this study was, therefore, to develop and demonstrate the DRAF, which automates these process steps as far as possible to make the methodology of mathematical optimization accessible to a broader user group. DRAF is meant to supply insights and decision support for academics in applied sciences, consulting engineers, and decision-makers in companies. Hence, it needs to be portable, easy-to-use, maintainable, editable, and extensible.

DR and Investments
The consideration of DR within optimal design planning generally results in higher capacities for assets and storage coupled to the electricity demand, since the associated lower average utilization rate is over-compensated for by the revenue or cost reductions from DR. Many analyses of DR potentials are limited to the potential of flexibility derived from existing assets. However, investment options that alter the existing flexibility, e.g., product storage extension for a flexible production process, or electrification measures, should also be considered when assessing the flexibility potential of L-MESs. E.g., Liu et al. [27] found synergistic effects when energy storage and DR of cooling, heating, and power were combined.

DR and Carbon Emissions
Electricity carbon emission factors (CEFs) can be categorized into grid-mix emission factors (XEFs) and marginal emission factors (MEFs). While XEFs are suitable for calculating carbon emission balances of energy consumption, MEFs are superior if the real carbon emission effects of a short-term demand change are to be approximated. Summerbell et al. [28] studied the cost and carbon reduction potentials of a cement plant through price-based DR using real-time prices (RTP). Although the carbon reduction potential was calculated from dynamic XEFs, they were not part of an objective function, so the carbon emissions could not be minimized. Baumgärtner et al. [29] calculated dynamic XEFs and MEFs for Germany and the year 2016 with an economic dispatch model and used them in the objective function of a multi-objective problem to design low-carbon L-MESs. They concluded that at the same costs, emissions can be reduced by 6% when using dynamic XEFs instead of annual XEFs and by up to 60% when using dynamic MEFs, which highlights the significance of using dynamic CEFs. In [30], optimization models considering DR have been reviewed. A lack of research on DR modeling for commercial and industrial consumers exists. Additionally, we strongly recommend the consideration of environmental effects by placing the carbon emissions in the objective function. In [31], dynamic XEFs and MEFs for 20 European countries were calculated, and a stylized price-based DR simulation was conducted. The results show that looking at national electricity systems, the effect of price-based DR on operational carbon emissions differs substantially from country to country and is dependent on the energy demand, generation mix, fuel costs, and carbon emission prices.

Multi-Objective Mixed Integer Linear Programming
We focus on multi-objective optimization, which combines two or more individual objective functions, e.g., the minimization of costs and carbon emissions, to gain a set of Pareto-optimal solutions [32]. These solutions are the best possible compromises and can be visualized as Pareto points on a scatter chart to give the user the option to choose from them.
The two most popular types of optimization methods are metaheuristics and mathematical programming [33]. Both types are used to optimize the operation and design of complex energy systems. Simplified, one can say that metaheuristics offer advantages for black-box models or for non-convex problems. In contrast to metaheuristics, mathematical programming can guarantee the achievement of a global optimal solution if an explicit equation-based model exists. Among the mathematical programming methods, mixedinteger linear programming (MILP) has shown to be effective for a wide range of cognate analyses. For instance, Zhang and Grossmann [34] listed 42 works with mathematical optimization models for industrial DR, among which, 35 were formulated as MILP models. For most energy-related components, the system behavior is nonlinear. However, the litera-ture largely approximates this system behavior via MILP models where binary variables are used within the piecewise linear approximations of nonlinear system dynamics [35]. Besides the approximation through formulating a MILP for a nonlinear system, there are further complexity reduction methods, such as time series aggregation, that can be applied to MILP models for energy system optimization [36].
In this work, we chose MILP, since the formulated optimization problems can be solved with global optimality within reasonable time frames. We assume that an explicit equationbased model exists and that all occurring non-linear phenomena can be approximately modeled as piecewise linear phenomena.

Open-Source
In the scientific context, the open-source idea is also becoming increasingly prevalent in the field of energy system analysis [37]. Publishing source code under an open-source license that is approved by the Open Source Initiative and listed under https://opensource. org/licenses (accessed on 29 June 2022) not only increases transparency and reproducibility, but also the quality of the software due to the increased incentive for collaboration [38]. In contrast to models intended for large geographic scales, models for potential analyses and planning of L-MESs are often created by consulting firms, where the release of source code is often not in line with current corporate practices [39]. As a consequence, tools for the analysis of L-MESs are less often published as open-source.

Other Energy System Frameworks/Models
In recent years, the scientific community has produced some comprehensive energy system modeling frameworks. In [40], 75 modeling tools were reviewed. A review of 24 energy system models and model generators was presented in [41]. Kriechbaum et al. [42] reviewed open-source modeling frameworks for grid-based multi-energy systems. A review on the concepts and validation models of multi-energy systems was conducted in [43]. Based on the reviews above, Table 1 presents a comparison between DRAF and other existing frameworks and model generators capable of building bottom-up models for operation and investment decision support of L-MESs. Note that the table contains some full-fledged frameworks, such as oemof, whose wide range of functions and aspects cannot be presented here. A powerful and versatile framework is the open energy modeling framework (oemof) [55]. More specifically, oemof is an organizational framework that bundles software for energy (system) modeling, such as the model generator oemof-solph [56]. However, oemof-solph uses Pyomo for the model generation, which builds models slower than the Gurobi Python interface GurobiPy [57]. Metadata such as units, parameter descriptions, and parameter value sources are not handled by the framework. While oemof cosmos offers several packages for parameter preparation, such as the load curve generation package demandlib [58], the preparation of important data for the modeling of price-based DR in L-MESs, such as dynamic CEFs, is not included. Very recently, Langiu et al., proposed the framework for component-oriented modeling and optimization for nonlinear design and operation (COMANDO) [48]. It focuses on nonlinear optimization; however, parameter preprocessing, interactive plotting, and metadata handling are not included in the package. There are also efforts to exploit the speed of the relatively new Julia programming language within energy system modeling, e.g., next energy modeling system for optimization (NEMO) [59], which was informed by the open source energy modelling system (OSeMOSYS). However, Julia is still a young programming language that lacks the massive community support needed for our target group, namely, programming beginners [60].
Some frameworks can be used in a flexible manner for manifold analyses (temporal and geographic resolutions) due to their modular structure and the clear separation between program logic and data, such as oemof [50]. This generic approach offers crucial advantages for large-scale power system analyses in terms of transparency, reusability, and maintainability. However, in order for the software to be used as decision support in an L-MES, it must be adapted to the specific problem with the following additional steps: • Model adaptation to industry-specific conditions; • Research and processing of market data, such as dynamic CEFs (depending on the country, year, and temporal resolution) and cost functions; • Generation of weather-dependent energy-relevant time series, such as energy yield time series for photovoltaics (PVs) or thermal load profiles; • Preparation, analysis, and plausibility checking of project-specific data, such as electrical load profiles, • Model parameterization; • Adaptation of result output functions, such as plots and tables to the particular data structure.
Commercial tools, such as the first three in Table 1, have addressed this need and can speed up the modeling process by specializing in specific application areas. A corresponding open-source alternative that leverages the potentials of the open-source idea (see Section 1.2.2) is currently not available.

Contributions
The main contributions of this paper are the development and demonstration of DRAF, an easy-to-use open-source Python decision support framework for optimizing DR-related design and operation of L-MESs. DRAF uses multi-objective MILP optimization to enable the user to quantify the cost and carbon emission reduction potential of existing and future flexibility options of L-MESs. The target groups of DRAF are researchers in the energy field and decision-makers in the commercial and industrial sectors. The source code of DRAF can be found at https://github.com/DrafProject/draf [61] (accessed on 28 June 2022). A secondary contribution of this study is the concise review of the current state of the art in open-source energy optimization software for DR.
The remainder of this paper is organized as follows. First, the most important elements and components of DRAF are described in Section 2, while referring to the extensive appendix with screenshots (Appendix A) and component templates (Appendix B). Second, the application of DRAF to three different simplified real-world case studies is demonstrated in Section 3. This is followed by a discussion, conclusion, and future research in Section 4.

The Demand Response Analysis Framework (DRAF)
2.1. Overview Figure 1 shows the main functional elements of DRAF. One can see that DRAF provides a toolbox for every step of typical energy system analysis and the optimization process of an L-MES decision maker. More specifically, DRAF is designed to answer the three questions illustrated in Figure 2.  Figure 1. Schematic depiction of DRAF's toolboxes, including examples and how they relate to the energy system analysis and optimization process and the used external resources/data. For example, parameterization is supported by DRAF through the parameter preparation tool TimeSeriesPrepper, which, e.g., provides electricity prices using the tool Elmada as an external resource. The general architecture of DRAF is presented in Figure 3. A typical use case is described in the following. First, a user imports and analyzes time series data of, e.g., electricity and natural gas historically purchased by the analyzed L-MES operator using Deman-dAnalyzer. Informed by the findings of the analyses, the user then instantiates a CaseStudy object of the desired analysis year and the country/address/coordinates of the L-MES. Subsequently, a first reference scenario that includes a model is added to the case study using component templates and the model generator. While DataBase is used here to provide and describe default parameters, the TimeSeriesPrepper prepares relevant time series, such as the day-ahead market prices, dynamic CEFs, or PV profiles. Different scenarios are then added by duplicating the reference scenario and changing specific parameters. After the model is solved by an external MILP solver, the results are stored in the cenario object. Finally, all parameters and results can be visualized either for each scenario (ScenarioPlotter) or for all scenarios in the case study (CaseStudyPlotter). The interconnected classes and modules allow for a fluent, explorative analysis process using the dot operator. E.g., cs. scens returns an overview of all defined scenarios; cs.scens.sc2.res.P_PV_fi_T.plot() plots the feed-in PV power of the scenario sc2. DRAF handles metadata, i.e., parameters can be stored together with descriptions, units, and sources. This motivates the input of metadata which can be used in plotting and exporting, prevents misunderstandings, and helps to document the meaning of an optimization model.  To carry out an analysis, the user initiates a CaseStudy and defines scenarios using the scenario generator, component templates, and parameter preparation tools (DataBase and TimeSeriesPrepper). The model generator builds optimization models from the model definition consisting of data and logic. After the optimization, the results can be plotted either via CaseStudyPlotter or via ScenarioPlotter using convenient dot notation (e.g., cs.plot.tables(); see Figure A3).
There are other energy modeling frameworks that also consider flexible demand. However, they do not focus on the modeling of L-MESs, such as individual commercial or industrial DERs. Since they are very generic, they technically allow the modeling of individual industrial sites; however, this abstraction comes with the cost of higher computational complexity and increased familiarization time with the tool for the user.
Additionally, they do not provide easy access to necessary market information, such as wholesale prices, dynamic emission factors, and technology investment costs. Versatility and portability are important aspects of DRAF. Both are ensured through a modular structure of DRAF that can be used as a generalized framework that can easily be adapted to specific applications.

Python as a High Level Programming Language
For the development of DRAF, the open-source, general-purpose programming language Python [62] 3.9 was chosen. A general-purpose programming language allows all steps of the analysis process from market data acquisition to data preparation, model building, scenario definition, and result visualization to be defined and reproducibly documented in a single environment. Since 2010, Python became the standard for new energy system models [41]. At the time of writing, Python is the most popular easy-to-use highlevel programming language [63]. It has rich libraries for data handling and visualization (e.g., Pandas [64], Matplotlib [65], and Plotly [66]). The open-source optimization modeling language Pyomo [67] provides access to important MILP solvers, such as Gurobi [68]. However, most importantly, Python has a vibrant community that supplies help and solutions to almost any problem that might occur.

Time Series Analysis Tools
As can be seen in Figure 1, the Time Series Analysis Tools refer to a DRAF toolbox. The included tools, DemandAnalyzer and PeakLoadAnalyzer, are described in the following.

DemandAnalyzer
Often, DR analyses are based on a time series of historic energy demands, e.g., the electricity demand of the year before. Before starting the modeling process, an analysis of this time series is helpful to validate the correct time-series length and to get key metrics, such as the peak-to-average ratio, load percentiles, and usage patterns. DRAF provides such an analysis with DemandAnalyzer. Figure A1 shows a screenshot of an example time series analysis.

PeakLoadAnalyzer
Based on the DemandAnalyzer object, the PeakLoadAnalyzer can be used; see screenshot in Figure A2. It shows the peak loads above a user-defined threshold and the cost reduction potential that originates from a given peak load price.

Parameter Preparation Tools
The Parameter Preparation Tools is a toolbox (cf. Figure 1) that supports the preparation of parameters for the optimization model.

TimeSeriesPrepper
DRAF's TimeSeriesPrepper allows the user to prepare time series such as the ambient air temperature; renewable electricity generation profiles; CEFs; and day-ahead market prices from basic input data, such as an address that converts to geographic coordinates, the analyzed year, and the time step width.

Carbon Emission Factors (CEFs) and Electricity Prices
In DRAF, CEFs and day-ahead prices for most European national electricity systems are automatically calculated for the given year and frequency with the open-source tool Elmada [69], as described in [31]. The latest Elmada version, v0.1.0 [70], supplies hourly and quarter-hourly time series for 30 European countries, mainly using data from the European Network of Transmission System Operators for Electricity (ENTSO-E) transparency platform [71]. For the CEFs, the user can choose between XEFs and MEFs, depending on the analysis question. For the electricity prices, the user can choose between day-ahead spot market prices which are referred to as RTP, time-of-use (TOU) pricing, and a flat price. In TOU, if not stated otherwise, time steps are grouped into high price and low price times; high price times apply between 8 a.m. to 8 p.m. during workdays. The high/low price is defined by the mean of RTP during high/low price times, respectively. The flat price is the annual mean of the RTPs.

Photovoltaic Power Profiles
If the user does not provide PV profiles for the analyzed location, DRAF uses the Global Solar Energy Estimator (GSEE) [72] to generate PV profiles from the geographic coordinates (or a valid address), global and diffuse radiation, and ambient temperature time series. For Germany only, DRAF checks for available weather data from the nearest weather station from [73] and downloads it in the background. This functionality may be extended to other countries as data become available.

Electrical and Thermal Load Profiles
The TimeSeriesPrepper module of DRAF also provides functions to create electrical and thermal load time series. In the electrical case, standard load profiles from [74] are used while considering public holidays from the given region with the Python holidays package [75]. For thermal load profiles, ambient temperature data from [73] are used to approximate heating and cooling demand time series.

DataBase
The reasoning behind the DataBase is to pragmatically provide the user with an expandable library of technical and market-related parameters, together with metadata such as units, descriptions, and scientific sources. Some of the values are used in the component template definitions in Appendix B.

Component-Based Model Generator
To avoid code repetition, and to make DRAF and the code written by the user maintainable, extensible, and adaptable, the toolbox Model Generator (cf. Figure 1) is implemented that creates model constructs in a lazy fashion; see Algorithm 1. Furthermore, the model generator keeps the model as light as possible, creating only the dimensions, parameters, optimization variables, and constraints that are needed. This allows the provision of adaptable and extensible component templates without limiting the user's freedom to build any other MILP model.
A component class consists of the two functions that define the component, param_func and model_func. The param_func defines dimensions, parameters, variables, and collectors for a given scenario. The model_func later uses these objects to build constraints and to connect the component to other components by contributing linear expressions to their collectors. The listings in Figure 4 show examples of these functions. The first listing defines a simple PV component. Note that dimensions and collectors are not needed for this simple PV component. The second listing shows relevant parts of the Main component, which defines general relationships that do not originate from a specific technical component.
For adding constraints and objective functions, the user can choose between Pyomo and GurobiPy syntax. While Pyomo supports different solvers, GurobiPy is limited to the commercial solver Gurobi but builds models with less computational effort.
Component interdependencies are considered using so-called collectors, which collect linear expressions across components and aggregate them in another component. C_inv_.values()). Collectors can collect scalar values, e.g, to aggregate investment costs of different components to total costs, or collect functions to access multi-dimensional vectors, e.g., to build an electricity balance for each time step; see Figure 5. If a component uses a collector, the constraints of that component must be built after the constraints of all components that contribute to that collector. This dependency restricts the order of submodel creation, which is resolved by executing a topological sort. This makes components reusable, so the user can conveniently choose from different storage and conversion technology components and modeling options, such as the consideration of investments or minimal part-load behavior, without inflating the model with overhead constructs. The user defines optimization models by using component templates (see Appendix B) and/or self-written technology components.

Component Templates
DRAF provides a set of component templates so that users do not need to start from scratch. The model is assembled from the Main component and multiple technology components. The Main component includes the objective function and all general sets, parameters, variables, and balances. Technology components are, e.g., energy demands, conversion and storage technologies, and interfaces to external entities, such as the electricity grid.

The Component Template Main
The Main component is a special component template that consists of the definition of the objective function and general balances and constraints. The Main component yields a deterministic multi-objective combined design and operation problem, which is described in the following. The formulation is partly based on previous work [76,77]. By default, DRAF works with equidistant time steps and assumes perfect foresight. The temporal resolution is defined by the user and not restricted; however, hourly or quarter-hourly resolution is required if the functionality of the TimeSeriesPrepper module is used; see Section 2.4.1. The default modeling horizon is one year, but can be customized. General sets are discrete time steps t ∈ T := {t 1 , . . . , t |T | }, components j ∈ J := {eDem, EG, PV, BES, . . . }, and flow types i ∈ I := {electricity, heat 1, heat 2, cool 1, cool 2, product 1, . . . }. In the follow-ing, optimization variables are denoted with bold symbols and are non-negative continuous variables unless stated differently.
The objective function of the MILP is to minimize the weighted sum of the total annualized cost (TAC) C tot and the annual carbon emissions CE tot : where α is the Pareto weighting factor ∈ [0..1] and π C , π CE are the Pareto normalization factors which can be identified by a simple algorithm to enable a more even distribution of the Pareto points; cf. Figure 14. X penalty j is a general penalty term which is only used in rare cases when there is an incentive that is not related to costs or emissions, e.g., to model uncontrolled battery electric vehicle (BEV) charging (in Appendix B.6), where the incentive is to charge based on time constraints.
The annualized total costs C tot are the sums of the annual operating costs, the annualized investment costs, and the maintenance costs of all components j ∈ J : The total carbon emissions CE tot are the sums of yearly operating carbon emissions of all components j ∈ J : The total annualized investment costs C invAnn j and the maintenance costs C rmi j are usually defined within storage and conversion components: C rmi j = k j,rmi c j,inv P j,capn where j stands for the component j, c j,inv are the specific investment costs, and P j,capn the new capacities. k af j (r, N j ) are the component-specific annuity factors defined in Equation (9) following, e.g., [78], where r is the calculated interest rate and N j is the operation life in years for component j.

Scenario Generation and Optimization
The scenario generator in DRAF provides convenient scenario generation. Scenarios can either be created manually or created in batches. For manual scenario creation, an existing object is cloned with sc = cs.add_scen(based_on=<scen_id>), whose parameters can be subsequently updated with sc.update_params(param1=value1, param2=value2 , ...); see also Algorithm 1. The batch scenario creation using cs.add_scens() can be seen as sensitivity analysis, which automatically creates a scenario for each combination of given parameters and parameter values. This is useful, e.g., for optimizing the system for different energy and carbon emission prices. When solving optimization models for a case study, the user can choose to solve the scenarios in parallel (cs.optimize(parallel=True)) using the distributed execution framework Ray [79] or serially to rely on the parallelization of the solver.

Visualization
DRAF provides a rich interactive visualization toolbox built into the CaseStudy and Scenario classes (cf. Figure 1). The dot notation allows convenient plotting since the data and metadata are internally fetched. E.g., after optimizing multiple scenarios, cs.plot. pareto() plots the Pareto front of all scenarios in the case study, similarly to Figure 14, and cs.scens.sc3.plot.sankey() plots the Sankey diagram of scenario sc3; see Figure A4. Interactive parameter and result exploration are available thanks to the diverse capabilities of Ipython [80] and Plotly [66]; see also Figure A5.

Case Studies
In the following, DRAF's features are demonstrated in three case studies that we consider to be of interest to the reader. The case studies are based on real companies in southern Germany with modified values for data protection reasons. In Case Study 1, the production schedule of a cement plant is optimized considering price-based DR. This case study was selected as it illustrates DRAF's support for flexible industrial production processes. In Case Study 2, the design of a battery energy storage (BES) and a PV system at an industrial site is optimized considering multiple flexibility applications and differentiating between existing and new technologies. Case Study 3 covers a more sophisticated superstructure for a greenfield L-MES. This last case study demonstrates multi-objective optimization, the value of the scenario generator, and the consideration of multiple temperature levels. The code for the presented case studies is available at https://github.com/DrafProject/draf_demo_case_studies (accessed on 28 June 2022). The calculations were performed using DRAF v.0.2.0 [61].

Case Study 1: Price-Based DR Potential of an Industrial Production Process
In the first case study, which is based on the previous work [81,82], we apply DRAF to the problem of quantifying the cost and carbon emission reduction potential of price-based DR of a cement milling process.
The setup of the case study is shown in Figure 6. In it, there are two electric cement mills that turn cement clinker into three different cement sorts which are stored in separate silos to serve a given cement demand. For each time step, the cement mills can either be powered down or produce one compatible cement sort in full-load with a sort and machine-specific production efficiency; see Figure 7 left. Machine 1 is compatible with sorts 1 and 3, and machine 2 is compatible with sorts 2 and 3; see Figure 6. The cement clinker supply is not a bottleneck in the production process, so it was modeled as unrestricted. The cement demand was generated by breaking down the total monthly demand into the working hours of the plant; see Figure 7 right. A MILP model was built using DRAF's component-templates Main, EG, pDEM, PP, and PS (see Appendices B.1,B.12-B.14) to minimize the TAC C tot of the system for 8760 hourly time steps from the year 2019. Therein, cement mills are represented by machines, cement clinker by raw material, and cement silos by product storage. Dynamic TOU and RTP pricing schemes and XEFs were prepared by DRAF's TimeSeriesPrepper, described in Section 2.4.1. The peak power price was e50 kW −1 p . A standard load profile for continuous production was scaled to the annual energy of 5 GWh and used as fixed electricity demand. The cement mills have nominal capacities of 3.5 MW each and a minimum part-load factor of 1. Each machine start-up and sort-change costs e10.00 due to the inefficient operation associated with it. Cement mill 1 was unavailable due to revisions from the 15th to 16th of March, and the cement mill 2 from the 15th to 16th of February. Each cement silo had a capacity and initial filling of 5 kt and a minimum filling level of 1 kt. Since part-load operation was not possible, the deviation between the last and the initial filling level was evaluated with a factor k and penalized via the operating costs. k was composed of the worst efficiency and three times the average electricity price of the year. Thus, the deviation was minimized without introducing infeasibility. For brevity, investment in storage extension was not allowed, even though this would be possible with the model and would be an interesting analysis.   Despite the technical constraints and the small price spreads of 2019 the load shifting led to savings of e149,000, (1.9% of electricity costs) and 1.4 kt CO2eq (9.2% of operating carbon emissions) per year. However, these are theoretical values, as complete foresight was assumed. The results for ten of the 365 optimized days are shown in Figure 8. The plots show that while the silo filling levels look similar, the scheduling of the cement sorts between the two scenarios differs substantially on an hourly basis. The electricity price is the main factor in the production decision. Exceptions are due to sort-switching and start-up costs. The Pearson correlation coefficient r between the electricity price and the purchased electricity for the whole year is −0.29 with a TOU pricing scheme and −0.58 with RTP; see Figure 9. Since machine 1 is more efficient than machine 2, the most energy-intensive sort 1 is only produced on machine 1. Sort 1 is only produced on machine 1 and sort 2 only on machine 2. The electrical peak demand was not lowered.   Figure 10 shows the setup and problem of Case Study 2.

Case Study 2: Design Optimization of a Multi-Use BES and PV System
For the inflexible electricity demand, anonymized real industrial data of the year 2020 were used, which are analyzed in Figure A1. Further input parameters were dynamic day-ahead market prices plus e62.3 MWh −1 electricity taxes and levies and e100 kW −1 peak electricity price. Specific investment cost forecasts for 2022 for PV and BES were taken from [83] with the values e384 kW −1 p and e209 kWh −1 , respectively. Figure A3 show the results of Case Study 2. Figure 11 shows the resulting electricity balance and the RTP of scenario optBesPv for one exemplary week (Monday-Sunday).  Figure 10. Setup and problem of Case Study 2. Setup: A company that can buy electricity from the grid and sell it to the grid has an existing 300 kW peak PV system, 1000 m 2 additionally available rooftop space for the installation of a new PV system, and an inflexible electricity demand. Problem: The design (nominal capacity/power) of the BES and the new PV is to be optimized assuming optimal charging and discharging of the BES considering peak shaving, RTPs through hourly wholesale prices, and the optimization of self-consumption.

Case Study 3: Multi-Objective Design and Operational Optimization of Thermal-Electric Sector Coupling
The third case study demonstrates the optimization of the design and operation of a more sophisticated industrial L-MES. This time it is a greenfield project-i.e., there are no existing technologies. Besides the inflexible electricity demand, the L-MES incorporates cooling and heating demands at different temperature levels. Figure 12 shows an overview and the superstructure of the analyzed L-MES. Up to 20 MWh/h of electricity can be bought from and sold to the grid with hourly prices and XEFs. Assuming plant-wide optimization with perfect foresight of XEFs, electricity prices, and PV yield profiles, Pareto-optimal design and operational configurations were to be identified that fulfilled the thermal and electrical energy demands. The following component templates from Appendix B were used: cDem, hDem, eDem, EG, Fuel, PV, BES, CHP, HOB, HP, P2H, H2H1, TES, and Main. The scheme in Figure 13 provides details on the modeling of the different temperature levels.

External energy supply
Electricity grid

Optimal control
Process heating demand 90°C / 70°C (H2)    Figure 13. Scheme with details on modeling different temperature levels for Case Study 3. The HPs could choose between three source and sink temperature levels. The evaporation and condensation temperatures were calculated assuming a 5 K temperature difference for heat exchange. The coefficient of performance was calculated from the evaporation and condensation temperatures. Assuming the installation of multiple HPs, multiple parallel operation modes, i.e., temperature combinations between evaporation and condensation, can exist; however, for the calculation of the annualized investment costs, the capacities of all HPs are aggregated, which significantly reduces the model's complexity. For more details on HP modeling, see Appendix B.10.

Storage levels
We modeled seven scenarios: The reference scenario REF and the six scenarios sc2-sc7. REF only allows heat-only-boilers (HOBs) and cooling machines that are modeled using the HP component by allowing only heat transfers from the cooling demands to the ambient temperature.
Scenarios sc2 to sc7 allow all technologies of the superstructure and all HP operating modes. They differ from each other only by the Pareto weighting factor α. The α values for the scenarios sc2 to sc7 were 0, 0.2, 0.4, 0.6, 0.8, and 1, respectively; i.e., sc2 optimized TAC (α = 0), scenarios sc3 to sc6 optimized Pareto efficiency (0 < α < 1), and sc7 optimized carbon emissions (α = 1). Figure 14 shows the TAC and annual carbon emissions of the resulting scenarios; and Figure 15 shows the resulting capacities, cost types, and distribution of the electricity exchange with the electricity grid. Due to the multi-objective optimization, scenarios sc3 to sc7 are Pareto-efficient-i.e., one objective value cannot be decreased without increasing the other. One can see that with increasing α values, i.e., increasing the focus on carbon emissions within the objective function, the investment cost and the annualized investment cost increase too. Compared to REF, sc2 and sc3 have lower TACs than REF by 36% and 9%, respectively, since higher annualized investment and maintenance costs are overcompensated by savings in operating costs. Additionally, sc2 and sc3 can reduce carbon emissions by 49% and 66%, respectively, compared to REF. Comparing sc2 and sc3 to REF produces no conflict of objectives, since REF is not on the Pareto frontier. The scenarios sc4 to sc7 have no economic advantage over REF; however, they do have an environmental advantage over REF. Scenarios sc3 and sc5 can be regarded as good trade-offs when looking at all available Pareto-efficient solutions. Scenario sc7 represents the highest possible carbon emission savings with 87%; however, the TACs are 41 times higher than in REF. Since in sc7 TACs are not considered, capacities were set to the highest possible values that are in place for each technology, e.g., 1 GWh for BES and 100 MWh for each TES. This is an unrealistic behavior that could be avoided by considering scope 3 emissions, which include the carbon emissions of the production of the energy technologies that are analogous to the annualized investment costs within the calculation of TACs. However, it demonstrates how scope 1 and scope 2 emissions can be reduced by increasing the system flexibility. As can be seen in Figure 15 bottom, in sc7 the technical upper limit of 20 MW of electrical power drawn from the grid is exploited during hours of low CEFs, which would stabilize the grid when there is a surplus of renewable energy. Screenshots of Sankey diagrams for REF and sc3 are shown in Figure A4.

Compared to REF:
TACs Emissions -87% Figure 14. Pareto-optimal configuration scenarios of Case Study 3. The dotted line approximates the Pareto frontier. REF is not on the Pareto front, since both objectives can be improved as, e.g., in sc2 or sc3. The broken y-axis was used to fit in the minimal-emission scenario sc7 (α = 1), which has a more than 41

Conclusions
We developed, described, and demonstrated DRAF, an open-source multi-objective decision support tool for L-MESs. By providing vital information about the environmental and economical potential of innovative energy technologies and services, DRAF lowers investment barriers of L-MES decision-makers. It considers load flexibility, energy efficiency, electrification, and their interdependencies in an integrated model without neglecting critical aspects, such as multiple temperature levels or DR of production processes. DRAF considers flexible electricity sources and sinks across the whole energy conversion chain of an L-MES, such as an industrial site. While providing useful pre-configured components that allow complex DR analyses with just a few lines of code, the tool setup does not restrict the users' freedom to build any possible MILP model. DRAF is needed since the existing software for potential identification of multi-energy systems is either (a) too generic to be practicably applied to L-MESs, (b) leaves out essential aspects such as temperature levels or the access to dynamic emission factors, or (c) is not open-source.
Three case studies demonstrate how different settings and applications can be modeled within DRAF. Case Study 1 shows how price-based DR of a production process can reduce costs and operating carbon emissions. Case Study 2 demonstrates a simple design and operational optimization problem for a PV and multi-use BES system. The more sophisticated design optimization of Case Study 3 demonstrates the consideration of multiple temperature levels, the selection of heat sources and sinks for the HP, and the results of a multi-objective Pareto analysis to select the optimal trade-off between economic considerations and the reduction of carbon emissions.
This paper shows only a small sample of the possibilities of DRAF. Future work is, therefore, the application of the framework for a detailed analysis and optimization of a real L-MES, such as an industrial company. The implementation of myopic and stochastic modeling in a rolling horizon fashion, tools for scenario generation and reduction, and the selection of typical days are also future work.

Appendix B. Component Templates Definitions
This section presents the mathematical formuation of component templates which can be classified as storages, conversion technologies, demands, interfaces, and a combination of them. For each technology, it contains an entity description table, the constraints, and the registrations to collectors. In the entity description table, all entities (parameters and variables) are listed with a source, unit, description, and default value(s) in the case of a scalar parameter. Thanks to consistent usage of naming conventions, the tables could programmatically be generated from the DRAF components, which also ensures consistency between the software and the paper. Furthermore, it demonstrates DRAF's possibilities in handling metadata such as units, data source information, and docstrings.
General parameters are: ∆ t is the width of the according time step in hours.  Φ source i=el,j=bes,t = P bes,out t Φ sink i=el,j=bes,t = P bes,in t C invAnn j=bes = 10 −3 E bes,capn c bes,inv k af (r, N bes ) C rmi j=bes = 10 −3 E bes,capn c bes,inv k bes,rmi Battery degradation is not considered. Please see [89] for more details regarding battery degradation in DR scenarios. A feed-in from the BES is not allowed. Following [90], η bes,ch and η bes,dis are calculated as the square root of the cycle efficiency, assuming them to be symmetrical.