pandapipes: an Open Source Piping Grid Calculation Package for the Application in Coupled Multi-Energy Grid Simulations

The increasing complexity of the design and operation evaluation process of multi-energy grids (MEGs) requires tools for the coupled simulation of power, gas and district heating grids. Most tools analyzed in this paper either do not allow coupling of infrastructures, simplify the grid model or are not publicly available. We introduce the open source piping grid simulation tool pandapipes that – in interaction with pandapower fulfills three crucial criteria: clear data structure, adaptable MEG model setup and performance. In an introduction to pandapipes we illustrate how it fulfills these criteria through its internal structure and demonstrate how it performs in comparison to STANET®. Then we show two case studies that have been performed with pandapipes already. The first case study demonstrates a peak shaving strategy as interaction of a local electricity and district heating grid in a small settlement. The second case study analyzes the potential of a power-to-gas device to serve as flexibility in a power grid under consideration of gas grid constraints. They both show the importance of a clear database, a simple simulation setup and good performance to set up different large and complex studies on grid infrastructure design and operation.


Introduction
In light of a new strategic view on the policy area "Clean Energy" within the European Green Deal [1] -a roadmap towards a sustainable economy developed by the European Commission -the energy infrastructure in Europe is expected to change drastically over the next years. On the one hand, the Hydrogen Strategy [2] is meant to promote the emergence of a hydrogen infrastructure that complements and partially replaces the existing natural gas infrastructure. The Energy System Integration Strategy [3] on the other hand aims at an increased interconnection of energy sectors and markets. It describes an energy infrastructure design that allows consumer needs to be satisfied through different energy carriers, thereby increasing its overall efficiency and resilience. These targeted changes will increase the complexity of the energy system as a whole.
As of today, different energy infrastructures are designed and operated independently from each other. Unlocking the potentials of sector integration requires a "new, holistic approach for both large-scale and local infrastructure planning" [3]. Methods and tools for infrastructure design need to reflect the increased complexity in order to address the many challenges on all scales and in all parts of the energy system design [4]. Challenges will occur in the following fields: greatly. For instance, Then et al. [5] give an overview of studies on the controversial role of the natural gas infrastructure and elaborate on the effects that a decline in energy demand can have on grid charges. In [6] they resume that increasing grid charges could accelerate gas grid defection as a self-induced effect. Kisse et al. [7] investigate the effects that changes in heating technologies for residential buildings have on investments into electrical and gas grid expansion and CO 2 -emissions by evaluating a number of scenarios. In order to improve the prediction of suitable supply scenarios for network operators in a multi-energy system (MES), methods and models have to consider all energy infrastructures in a combined optimization approach.
• Grid integration and expansion in coupled energy infrastructures: The progressing electrification and renewable generation increases both load and generation in distribution grids. All new devices and power plants have to be integrated on the respective voltage level which often urges network operators to expand or re-design their power grids. This planning process in future needs to address the increasing interconnection to other infrastructures and new opportunities for flexibility trading or storing energy [8]. Approaches should take into account the possibility to rededicate parts of the infrastructure to other carriers or even remove infrastructure that is not required anymore [2,5]. It is crucial to identify suitable investment paths towards a future supply infrastructure taking into account the longevity and high capital expenditure of infrastructure assets. These challenges also have to be reflected by large scale grid integration and expansion studies like [9,10]. • Analysis and optimization of operation strategies: Sector coupling facilities are expected to deliver ancillary services for power system operation (e.g. frequency support, voltage support, congestion management etc.). Although such services are mainly required in the power grid, operation strategies should be optimized by respecting the operational constraints of all connected infrastructures. For example, Liu et al. [11] analyze different operation modes of combined heat and power (CHP) plants and check which one is most suitable to serve the power grid without violating constraints in the gas grid. In [12] an optimal power flow (OPF) implementation is presented that includes constraints from the gas grid as side conditions. • District planning: The planning phase of new districts offers great potential for low-cost decarbonized energy supply by implementing small MESs with optimized grid infrastructure [13]. In such small systems, choosing and optimizing the operation strategy is crucial in order to use energy locally or to offer ancillary services to the external power grid [14]. • Market design: Sector coupling facilities are well suited to deliver flexibilities and other ancillary services, for which they also need to be remunerated. In [15] Mancarella et al. suggest an approach to determine the profitability of MES to deliver ancillary services by considering the revenue and the cost of energy shifting. With an increased interest in local energy markets [16] comes a need to analyze approaches with respect to the technical and economic performance using models that comprise the complex state of the whole MES and market mechanisms.
The challenges described above require powerful tools for coupled MEG simulation and optimization. In this paper we introduce the open source tool pandapipes, which is meant to simulate stationary and quasi-stationary flows in pipes, considering hydraulic and thermal properties of the flowing fluid. It complements the open source power system analysis tool pandapower [17], as it is based on the same structure and implements similar core methods, such as the Newton-Raphson solver. The two tools can be combined into a powerful and easy-to-use MEG simulation environment with the help of the pandapipes multi-energy module. As both tools are published under the terms of a 3-clause BSD license, they can be used by anyone for any purpose free of charge. Being based on python, these tools can be combined with any other python package for further analysis. For instance, there are interfaces to NetworkX [18] for topology analysis, to matplotlib [19] for plotting and to geopandas [20] for processing geographical information. All these characteristics enhance the potential of pandapipes to automate processes and easily adapt the required models and simulation setups. Therefore pandapipes is well suited to address the challenges of energy infrastructure analysis and planning. Since the coupling with pandapower is very simple, many use cases that have already been addressed with pandapower can be adapted to an integrated approach.
In the present paper we want to illustrate the innovative potential of pandapipes. Section 2 gives an overview over different approaches to MEG simulation and derives three criteria to be addressed by a dedicated tool: clear data structure, adaptable MEG model setup and performance. Section 3 introduces the tool and describes its structure and model implementations, thereby illustrating how these criteria are addressed. Section 4 demonstrates two use cases that were implemented with the help of pandapipes already. In the first use case, we analyze local peak shaving strategies with the help of a heat pump connected to a district heating grid. It is shown that by applying the introduced combination of pandapower and pandapipes, operation strategies and their effect on all involved networks can be examined. In the second use case, we analyze the potential of power-to-gas (P2G) devices to reduce congestions in a power grid without violating constraints of the connected gas grid. The presented simulation approach enables large scale analyses on the effectiveness of P2G devices in different configurations to find a suitable system design. Section 5 concludes on how pandapipes enhances the toolbox of researchers working in the field of MEG analysis and gives insights into possible enhancements to further simplify MEG simulations and to address other types of problems.

Multi-Energy Grid Simulation: Overview of Approaches and Requirements
Approaches for the simulation of coupled energy infrastructures can have different levels of detail in specific areas, making them especially suitable for specific use cases. Thus a large tool landscape has evolved of which Table 1 shows an excerpt. We characterize the tools by licensing, modeled infrastructures, model detail and the possibility to conduct coupled simulations. Of special interest for the research community are open source tools, as they are customizable, automatable and mostly free of charge. Commercial tools usually have a strong focus on the user interface and usability. We distinguish three tool categories: grid simulation tools, energy system modeling and optimization tools and multi-energy grid (MEG) simulation tools.
The characteristic of grid simulation (GS) tools is that they offer detailed grid models, sometimes even for different infrastructures, but they cannot be coupled within the simulation. Energy system modeling and optimization (ESMO) tools analyze power flows between defined areas with a strong focus on energy balancing. Therefore they precisely model power plants and devices with their control strategies and characteristics like efficiency or ramp rates. As they do not provide detailed grid models, they can easily integrate all infrastructures within one simulation. Such tools can be used for optimal power system design or power plant deployment planning, i.e. in a flexibility market model used for power balancing.
None of those tools offers sophisticated models for integrated infrastructures, as either the grid model is simplified or a coupled simulation is not possible. For the simulation of MEGs, two approaches are prominent: Combining dedicated tools with the help of a co-simulation platform and integrating all models in one tool. As pointed out in [21], the main advantage of a co-simulation approach is the possibility of integrating any tool while leaving the development to the respective experts. Furthermore distributed simulations can be performed online so that participants can exclude exchange of internal data. But there are also downsides to this approach: Firstly exchanging data between two or more tools via an external interface is time consuming and usually acts as a bottleneck for complex simulations resolved in space and time. Moreover, the design, functionality and licensing of coupled tools might introduce difficulties for the user to set up the model. Most co-simulation approaches are tailor-made [22] and those with a wider scope are still not fitting for specific applications [23]. In the last years, also dedicated tools were developed to address problems in the field of MEG simulation and to overcome the above stated problems of co-simulation. Such tools that enable the user to easily run simulations with models for all infrastructures are rare and have not yet found their way into the open source community, as can be seen from Table 1.
The main challenge of setting up MEG simulations lies in the complexity of the addressed problems. Therefore a tool that is designed to solve a multitude of different problems in this field has to focus on tackling this challenge, especially considering that users' experience might differ greatly. A consistent design is crucial and the first simple models should be very easy to build, while complex models should be solved efficiently. From these requirements we derive three criteria that should be addressed by a dedicated tool in our opinion: • Clear Data Structure: Mastering the complexity of MEG simulations is only possible, if the data is contained in a clearly structured database that allows storing large amounts of data. This is true for the model data of the coupled grids and the output data, especially from time series simulations. When performing large scale analysis studies, the simulation time plays an important role which has to be considered in the design of the tool. But even for daily routines spending a lot of time waiting for calculations to finish can be very inconvenient and inefficient.
As most of the MEG simulation tools introduced in Table 1 are closed-source or based on proprietary software, these criteria are hard to analyze for them, but some publications give hints. SAInt implements an OPF formulation considering AC power flow equations, transient gas grid equations and coupling terms which can be used to analyze outages in coupled transmission grids. A lot of emphasis is put on an increased performance by using a specialized solver. In HyFlow active power flow calculations are performed with models for power transport gas and heat grids, linked by simplified coupling terms and solved with a Newton-Raphson approach. The field of application is the analysis of power balancing in a cellular approach. MYNTS uses an expression tree formulation of component equations for automatic differentiation, which are solved with non-linear programming solvers. Its main field of application is the evaluation of supply scenarios for gas transport grids, but also coupling terms can be implemented with the generic approach. TransiEnt is a modelica library that defines transient and steady-state equations for the calculation of MES. Its modular approach and different levels of simulation detail allow different simulation setups. Applications have shown the interconnection of operation strategies for different components.
Next to TransiEnt, pandapipes with its specialized multi-energy module is the only MEG simulation tool with an open-source license. Unlike the former, it can even be used without any proprietary software required. It was developed along practical use cases and designed to fulfill the three aforementioned criteria. The multi-energy module integrates pandapipes and pandapower to address challenges in the field of coupled MEG planning and operation. pandapower has already shown to be well suited as a tool for solving automated planning problems [40], thereby considering different effects like restoration [41], operation strategies [42] or asset management [43]. It was also used for testing new approaches to contingency analysis [44], curtailment minimization of distributed generators [45], state estimation [46] and fault occurrence evaluation [47]. Its special potential lies in the application within studies that automatedly analyze large grid areas, such as grid integration studies [10,48] or energy loss studies [49,50]. Due to its similar design and performance, pandapipes is capable of performing similar studies with respect to piping grids, and in combination with pandapower also to MEGs.

Overview of pandapipes
In the previous section, we introduced three criteria for a MEG simulation tool. In the following we illustrate how pandapipes meets those criteria through a comprehensible design and an efficient calculation method. We introduce the pandapipes grid data structure and the controller architecture, then demonstrate how this architecture can be used to build up multi-energy grid simulation models and finally compare the performance of the pipeflow calculation with calculation times of the popular proprietary software STANET®.

Illustrating Clear Data Structure: Architecture of pandapipes
In pandapipes, nodes and edges of a network are defined by component models. The most important components are junctions, defining nodes, and pipes, which define a connection between two junctions. Every component model introduces equations, describing the physical properties that are solved for. The equations are assembled into a nonlinear system of equations, solved with the Newton-Raphson method. As a result, pandapipes provides pressure values at all junctions in the network and the corresponding flow velocities along the different pipes. If a heat grid is modeled, pandapipes also determines temperature values at the junctions. Further properties can be derived from these dependent variables.

The pandapipes grid structure
There are different parameters which have to be provided by the user to properly set up the mentioned system of equations. This data is stored in a dedicated data structure, shown in Figure 1. The pandapipes net is a container for different simulation input and output. It contains pandas tables for different components and additional parameters required to run a simulation, which is an analogy of the pandapower structure. Utilizing the pandas format allows a simple creation of components, such as pipes, valves or pumps with nameplate parameters. There are two main types of tables: Component tables contain all the modeled grid elements with their respective nameplate parameters. For example, the pipe table contains the two connected junctions along with the parameters length, diameter, roughness, an additional loss coefficient and a flag for its activeness. Those values are required for the hydraulic state simulation. Furthermore it is possible to set an external heat source or sink, an external temperature and a heat transfer coefficient for the heat transfer calculation. Result tables contain the simulation results for the given component. For some components it is also useful to record geo-information which is then stored in another table.
Properties of the operating fluid, such as density for different temperatures and pressures, can be freely defined by the user or chosen from an included fluid library and will be stored in the grid container as well. Other stored parameters are standard types for components as pumps or pipes and further options to be set by the user. The available component models are held in a list for each pandapipes net for the internal process. Other internal data include the internal array structure and lookups for the conversion between the component tables and the internal structure.

The pipeflow procedure
After defining the input data inside the component tables, the calculation can be started by calling the pipeflow function. Figure 2 shows the procedure of the calculation. With the help of the component models, all component tables are transferred to the internal array structure and equations are introduced for the derivatives of the state variables. A connectivity check makes sure that disconnected grid areas are excluded from the calculation. Afterwards the hydraulic state variables are solved for with the help of a Newton-Raphson solver which usually requires several iterations to converge. For heat grids, a subsequent heat transfer calculation also uses a Newton-Raphson solver to calculate the node temperatures. In the end, the state variables and derived results are transferred to the result tables.   The most important equations are introduced by branch components which define the pressure loss and the heat transfer model. These models differ between compressible and incompressible media. In the current version of pandapipes it is possible to calculate hydraulic properties for compressible and incompressible media and to perform a subsequent heat transfer calculation for incompressible media. An overview of the most important equations is presented in the Appendices A and B.

The controller architecture and time series simulations
Typically not only one stationary state is of interest, but instead the grid state has to be observed over a specified period of time. For this purpose pandapipes provides a dedicated function that starts a loop over the time period, as shown in Figure 3. Connected devices can be modeled with the help of time series data stored in a file or through a dedicated control scheme. Time series data are loaded in every time step and handed over to the grid structure in the time step initialization. For the implementation of control schemes, a controller architecture is implemented in pandapipes in analogy to pandapower (cf. [51]). Controllers can be used to control process variables (e.g. the pressure at a specific node) by setting reference variables (e.g. the feed-in at this node). As they are implemented as python classes, a new controller class can inherit properties, especially the interface to the data structure, from a controller parent class. Developers who need to introduce a new model can thus focus on the physical modelling. In every time step of the time series loop, an internal control loop is started that iterates over the controllers and performs several pipeflows until all controllers are converged.

Illustrating Adaptable MEG Model Setup: Introduction to pandapipes multi-energy
The coupling points of energy infrastructure are facilities like CHP plants, heat pumps or P2G devices. As their specific models are usually not inherent to a grid calculation, they need to be implemented in a separate structure. In literature the concept of energy hubs is prominent to model such facilities [52,53]. An energy hub is a generic formulation of a unit that stores or converts energy between different infrastructures. The energy conversion and storage must follow rules which are usually defined by a control scheme.  For this purpose the multi-energy module of pandapipes extends the controller architecture to couple different grids with the help of special coupling controllers. It defines a superordinate MEG structure which contains several grids and multi-energy controllers that can model the energy transfer between them, as shown in Figure 4. For example, a heat pump can be modeled as a unit that converts electrical energy as an input from a power grid into a heat flow as an output to a heat grid. With the help of a heat pump characteristic and the heat grid temperature as second input, an operating point with its coefficient of performance (COP) can be calculated. If it is used to control the temperature at the outlet node, several control iterations with successive pipeflows can be necessary to set the power correctly. Possible models with the exchanged data of a heat pump and a P2G device are shown in Figure A1 of Appendix C. As the controllers exchange data with the connected grids for the purpose of controlling certain variables, it is easily possible for them to also ensure a correct model of the energy transfer. This way they serve as energy hub model implementations as well.
The described approach differs from those used in other MEG simulation tools, as the equations of components in different infrastructures are not summed up in a single system matrix or problem formulation. Therefore it allows for a high degree of flexibility; however it leaves power balance checks to the controller implementations and might lead to a slightly lower performance. Although different calculations for different infrastructures are necessary like in an approach using a co-simulation platform, all calculations are performed within the same environment and no communication interface between different tools is required. This drastically reduces the communication overhead and simplifies the model setup. Moreover, the co-simulation is only performed on the inner structure of the model; a user does not need to learn the specifics of different tools that are plugged together, as pandapipes and pandapower implement an analogous user interface.
In the case of MEG simulations time series studies are of special interest, as the coupling facilities can be used to shift power peaks in time. Time series simulations also allow to analyze the effectiveness of a certain control scheme. For this purpose, the multi-energy module defines a specialized time series simulation setup for MEGs. In each time step the controllers of the individual grids are called along with the multi-energy controllers for coupling. Several simulations of each grid might be  required until all controllers converge. The multi-energy grid model is the core of the MEG simulation environment based on pandapipes and pandapower. This simulation environment can serve as basis for a framework to tackle challenges in the field of planning and operation analysis of MEGs, as shown in the right part of Figure 4. Dedicated modules might address the fields of operation optimization, automated planning and topology optimization, studies on ancillary services and an integrated power flow optimization.

Illustrating Performance: Comparison between pandapipes and STANET®
Studies that analyze the operation and design of MEGs usually require a large number of evaluations in different configurations. Furthermore, these studies often require time series simulations as mentioned previously. Therefore, a detailed study of different configurations of a MEG setup can easily require several thousands to millions of single simulation steps, demanding a highly performant core. This requirement can be satisfied by the use of pandapipes. In order to analyze the performance of pandapipes, we compare it to STANET®, one of the leading piping grid simulation tools available on the German market.

Simulation Setup for the Performance Comparison between pandapipes and STANET®
As performance comparisons are highly dependent on the given setup our goal is to create an environment with an equal base for both tools. The following calculations are all run on an ordinary laptop with the specifications given in Table 2. For the performance comparison, we conduct time series calculations for one entire day with a resolution of 15 minutes (96 time steps) for three grids of different size. We repeat the calculation 20 times to minimize the effect of outlier results. The following three grids, shown in Appendix D, are analyzed: • The district heating grid ( Figure A2), consists of four heat exchangers. The demand of each heat exchanger is constant over time. No fluid supply is required, as the net is considered to be a closed system. Solely the pump ensures a circulation of the fluid and provides the required heat supply.   The intention of creating three different, quite diverse grids is to cover many aspects affecting the performance in an efficient way, however, being aware that this is by no means exhaustive. One of these relevant aspects is the considered energy carrier. For example, simulating a district heating grid in pandapipes always requires a subsequent heat transfer simulation in addition to the hydraulic simulation, as can be seen in the flowchart in Figure 2. This is usually not necessary for simulations of water and gas grids. Therefore, all three simulation options which can be modelled with pandapipes are represented in this paper: hydraulic calculation of compressible and incompressible media and the additional heat transfer calculation of incompressible media. Another aspect we consider as relevant is the number of different components integrated in a grid model. Therefore, the water grid comprises all components that can be modelled in pandapipes, except for the heat exchanger. Also the size of the grid has a significant effect on the simulation speed. Therefore, we steadily increased the size from the district heating grid up to the gas grid. A last effect we identified is the meshing degree. Therefore, the district heating and water grid have up to 30 % higher meshing degrees compared to the gas grid.
An overview of relevant grid characteristics is given in Table 3.  Table 3 shows that the water grid has the highest number of different component types. The district heating and gas grid contain the same number of component types, yet the types differ as well as the grid size. The meshing degree of the water grid is slightly higher compared to the district heating grid and much higher compared to the gas grid.
The table also shows that the total number of components installed in the different grids increases at a disproportionately low, whereas the number of profiles increases at a disproportionately high rate with grid size.

Comparison of Calculation Time
The performance results are displayed in Figure 5, retrieved from 20 simulation runs for each grid conducted in pandapipes and STANET®, respectively. The average results are represented by the colored bars, while the standard deviation is visualized by a black line.  We distinguish between the pure simulation time and the total time. Whereas the simulation time just comprises the required time for the pipe flow, the total time additionally includes the overhead time like loading profiles and saving results. In each bar plot the average simulation time is represented by the lower number, the total time by the number above.
Both tools have in common, that a grid size increase results in a rise of the total time both for pandapipes and STANET®, however, with different characteristics.
In case of STANET®, both overhead time and pure simulation time increase with size. pandapipes, in contrast, reveals a similar overhead time in case of the heat and water grid, while depicting a sudden jump during the gas grid simulation. Also the pure simulation time shows a peculiarity, as the solver requires more time for the smaller water grid than for the much bigger gas grid. This anomaly, however, can be led back to the number of solver iterations to reach convergence. While in case of the gas grid, usually two to three iterations are required, the solver is called seven to eight times in case of the water grid.
A comparison of both tools reveals that pandapipes is always faster than STANET®. While in case of the smallest grid, pandapipes is only 3 times faster, the difference becomes more predominant with grid size making pandapipes up to nine times faster in case of the gas grid. One reason can be found in the additional overhead like the graphical user interface (GUI), which STANET® provides.
However, comparing the pure simulation times with one another, the solver alone also claims more time in case of STANET® compared to pandapipes. One reason, we identified, might be the higher accuracy in case of STANET®, what also can be seen in the absolute deviations between the pandapipes and STANET® results, visualized in Figure 6. While STANET® for example considers the temperature dependency of the dynamic viscosity, we accept a slightly smaller accuracy for the gain of a much better performance. Another reason might be the readout speed of the profile data. In our examinations, the time spent to read the data from the disc is almost negligible. This matter of fact, however, shifts as soon as another data format like .csv is used, slowing down pandapipes massively.  Figure 6. Absolute deviations of the pandapipes compared to the STANET® results. From left to right, one can see the absolute deviations of pressure at the nodes, fluid velocity in the pipes and temperature at the junctions for each grid, respectively. As only the heat grid considers a heat transfer calculation, we solely displayed for this grid the absolute temperature deviation. In order to classify the deviation we also calculated the mean pressure and the mean absolute deviation and set them in relation. The overview can be seen in Table 4 Preprints (www.preprints.org) | NOT PEER-REVIEWED | Posted: 6 October 2020 doi:10.20944/preprints202010.0132.v1 Also other reasons can be causal for the diverging simulation speed including system dependencies or the investigated grid, both we only analyzed to a limited degree. Furthermore, our look below the hood of STANET® was only possible as far as the manual took us, as the software itself is closed source. Therefore, our conclusions are limited by the information in the instruction book. However, in all our investigations we tried to be as transparent and unbiased as possible. In the Innonex project [56], operation strategies were developed to convert and store available excess power generated by photovoltaic (PV) plants as thermal energy. The excess energy supplied a centrally aligned heat pump connected to a district heating grid, as well as electrical heaters installed in the drinking hot water storages of the consumer households. In this way, the excess power could be used to cover peak loads occurring at a later time. A detailed model of the district heating grid made it possible to evaluate temperature levels and thus thermal losses in the grid. However, no model was created for the electrical grid, which made it impossible to react on certain events like line overloadings. Thus, the state of the electrical network was not respected by the operation strategy.

Solving Problems of Coupled Multi-Energy Grid
A similar, but more complex simulation setup consisting of different tools combined in a co-simulation framework and coupled with an optimization tool was used in [57] and [21]. The tools are used to maximize self-consumption from RES and minimize CO 2 -emissions and electricity imports to the district. Typically, a co-simulation approach acts as a performance bottleneck. This is different with pandapipes multi-energy where all models are included in one tool, thereby simplifying the model setup.
In the use case that we analyze, the heat energy -consisting of the demand for space heating and drinking hot water -for a small settlement is supplied by a district heating grid fed by a centrally positioned heat pump. The operation of the heat pump depends on the status of the power grid. Therefore, a coupling between the electrical and district heating grid has to be provided. It is assumed that the heat pump has access to a low temperature heat source, so that an inlet temperature of 45°C or alternatively 55°C can be supplied. However, regulations demand a minimum temperature of the drinking hot water of 55°C.
Every building connected to the district heating grid has a drinking hot water storage. The storage extracts heat from the heating grid using a heat exchanger. To lift the temperature from the district Preprints (www.preprints.org) | NOT PEER-REVIEWED | Posted: 6 October 2020 doi:10.20944/preprints202010.0132.v1 heating grid temperature level to the required temperature of at least 55°C, an electric heater is placed inside the storage. All buildings are equipped with PV plants whose primary purpose is the electrical supply of household appliances. If excess power is available, it can be used by the electric heaters inside the storages. Further excess power is then fed into the electrical distribution grid.
The heat pump heats its storage to 55°C if the prosumer households provide sufficient renewable power to cover the heat pumps electrical demand. Otherwise the provided heat pump temperature is 45°C and the required power is supplied partially by excess PV feed-in and power drawn from the connected distribution grid. As an additional constraint, the temperature of the heat pump can only be raised if this operation mode does not lead to a line overloading in the connected electrical grid. If a line overloading occurs even at a provided temperature of 45°C, the heat pump enters a power saving mode, where nearly no mass flow is provided and no power is consumed.
The settlement consists of 12 buildings. The network topologies for the power and district heating grid are shown in Figure 7. A detailed flow chart of the implemented control strategies for the electrical heaters and the heat pump can be found in Appendix E.

Use Case Implementation
Both pandapower and pandapipes only have access to a small set of components like sinks and sources. More detailed components, like the required prosumers and the heat pump, have to be modeled as controllers. Figure 8 shows a sketch of the prosumer controller and its components. As depicted, each water storage is connected to the pandapipes district heating grid via a heat exchanger. In every time step, the mass flow and inlet temperature (T net ,ṁ) are extracted and used as an input to simulate the water storage temperature. If activated, the electric heater is powered with a constant value of 1 kW.  In addition, the amount of extracted heat by the water storage heat exchanger and the heat required for space heating is transferred to the pandapipes net (q net ).
Besides, the controller is connected to a sgen component in the power grid model. Depending on the sign of the power variable, this component may represent a producer or a consumer. In each time iteration, the current magnitude of excess power is determined by the prosumer controller and transferred to the sgen component.
The available excess power is computed according to equation 1, where PV(t) corresponds to the power provided by the PV plant and load(t) corresponds to the electrical demand of all devices in the household except the electrical heater heating the water storage. The latter is described with the variable dhw el (t). If ep(t) is negative, the available PV power is not sufficient to cover the electrical demand. In this case, power has to be fed in by the power grid.
The described data exchange between the controller and the connected grid models is started via the controller's time step initialization, where models of arbitrary complexity may be included. The water storage temperature T is computed via the differential equation 2 A numerical solution for the next storage temperature is determined using a differential equation solver provided by the python package scipy.
The time series for the PV power plant, the electrical household loads and for space heating are provided as external sources. They have been generated with models described in more detail in [58] and [59].
A diagram of the implemented heat pump controller and its interfaces to the pandapower and pandapipes components is shown on Figure 9. It also has access to a water storage. The heat pump feeds directly into the storage (with a provided temperature T in ), which in turn is connected to the pandapipes district heating grid. The current storage temperature (T storage ) computed in every time step as well as the mass flow (ṁ) provided by the pump are parameters for the pandapipes component. The heat pump is also connected to a load component of the power grid. The controller transfers the required power (P el ) for operating the heat pump in every time step.

Result Evaluation
We examine two days in July of the selected test reference year 2015 [60]. In this month, the excess power produced by the PV power plants is typically enough to supply not only the electrical heaters of the buildings but also the heat pump. Since space heating is hardly required, it is not discussed in the following. On the first day, the blue power line depicted on Figure 7 would have been overloaded, so the storage temperature does not rise to 55°C. This is only possible on day 2. Figure 10 does not only show the temperature of the heat pump storage, but also a plot for the heat pump control signal. Additionally, the excess power is displayed.  The blue curve shows the available excess power after the household loads of each house has been subtracted. The orange line represents the total excess power after the power requirements for the electrical heaters are subtracted. The green line shows the amount of total available excess power. We see that the amount of available excess power is further reduced by the heat pump operation. Remaining excess power could be fed into the external grid.
A similar result focused on one building in the settlement is shown in Figure 11. The upper panel shows the storage temperature during the two days. The required thermal power is supplied by the district heating grid and the electric heater. The slope of the storage temperature is connected to the demand profile shown in the plot in the middle. The lower panel displays the available excess power for the observed household. Again, the blue curve represents the generated PV output minus the household loads and the orange curve adds the required heater power. Note that the storage temperature plot reflects the control algorithm of Figure A5. For example, the heater is turned off after the secondary setpoint temperature of 70°C is reached.

State of the Art Review
In some regions of Germany the fluctuating feed-in from renewable energy sources (RES) already leads to congestions on the distribution grid level. One possibility to mitigate such congestion situations is to expand and reinforce the distribution grid, which is a long-term strategy with high capital expenses. Other options include curtailment of renewable feed-in, redispatch with other power plants or demand-side management. Those operational interventions are short-term and the single measures are rather cheap. In order to implement them in a cost-efficient way, market platforms for trading flexibility options get more attention recently [61]. sector coupling facilities like CHP plants or power-to-gas devices are considered to have great potential as such flexibility providers [62,63]. Preprints (www.preprints.org) | NOT PEER-REVIEWED | Posted: 6 October 2020 doi:10.20944/preprints202010.0132.v1 Since local flexibility markets are a rather new concept, the model development is still in its early phase. In [64] CHP plants are modeled with a detailed bidding strategy and some operational limits, but the coupled gas grid and heating infrastructure are neglected completely. In [65] an integrated gas and power grid optimal power flow is presented to use CHP plants with thermal storages as flexibilities for operation cost minimization. The authors show its performance in a small test case for a two-day period. An analysis of linepack flexibility is presented in [66] with the help of a multi-stage procedure including several separate evaluations of the gas and power grid. The potential of P2G devices to reduce necessary curtailment of RES is analyzed in [67] with the help of an integrated optimization formulation for the power and gas grid, modeling P2G devices with a constant efficiency.

Use Case Implementation
A model for analyzing flexibility provision was also implemented in pandapower within the project NEW 4.0 [68]. This model is similar to the approach in [69] and will be described in detail in separate publications. It can be used to find optimal operational strategies for flexibilities and analyze the relationship between investment and operational measures for the distribution system operator. An OPF is used to set the active and reactive power values of flexibilities to a cost optimal value which deviates as little as possible from their predefined set points. For P2G devices, a feedback from the coupled grid should be taken into account. There the P2G devices' operation could lead to the violation of operational constraints, thus decreasing their potential to serve as flexibilities in the power grid.
An approach to integrate gas grid constraints into the flexibility provision model was developed with pandapipes multi-energy on the basis of synthetically created grids. The used medium voltage electric grid [51,70] geographically overlaps with the medium pressure gas grid with direct household supply [7], as shown in Figures 12 and 13. In order to better show the effects occurring in the gas grid, we altered the gas grid data by shifting the slack node, lowering the overall pressure level and artificially increasing the pipe resistance. The grid data is complemented with time series for load and generation. Standard load profiles are used for electric loads and load profiles based on typical heating demand considering probabilistic occupancy are used for the gas consumers [71,72]. Time series of PV and wind power feed-in are provided for all RES sites and combined to three scenarios: the PV scenario with 100% PV feed-in, the wind scenario with 100% feed-in from wind power and the mix scenario with 50% feed-in from PV and 50% from wind power. Table 5 gives a small summary of the maximum power and annually generated or used energy.
The operational constraints that lead to adaptations of the RES and P2G feed-in are shown in Table 6, along with boundary conditions. We consider a maximum line current and node voltage in the power grid and a maximum pressure and P2G feed-in in the gas grid. As we analyze two different P2G processes -one with successive methanation and one with direct hydrogren feed-in -there are two different constraints that limit the feed-in. In the configurations with methane feed-in we assume that there is a uni-directional flow at the slack node, because normally gas pressure regulation stations do not contain a compressor that could feed in gas to the high pressure grid. In the configurations with hydrogen feed-in, we assume that the volume fraction at the connected node must remain below 10%.
We perform a simulation for a full year in 15 minute resolution with the approach depicted in Figure 14. This approach considers differences in demand and weather occurring over the course of the year. It makes use of the special structure and the time series and controller architecture defined in pandapipes multi-energy. In each time step constant controllers write values from load and generation profiles to the respective pandapipes and pandapower grid tables. Based on the results of a subsequent power flow, a central power grid controller checks for constraint violations according to Table 6. In case of violations, an OPF can reset the power of the P2G devices and RES to comply with the power grid constraints. The next triggered controller is the P2G controller. It first transforms the adapted power of the P2G devices to a mass flow according to a P2G model and performs a pipe flow simulation.    based on the total gas consumption and the necessary mass flow reduction is distributed over all P2G devices proportionally to their identified feed-in. If the maximum pressure constraint is violated, the P2G devices are replaced with slack nodes with fixed pressure at the maximum value to identify the allowed feed-in. After transforming the maximum mass flows back to the maximum powers of the P2G devices, they will be considered by the OPF in the next control loop. This feedback from one controller to the other is not a realistic control loop and thus can only be used to identify the technical potential. In reality these calculations have to be performed with forecasted grid states introducing an uncertainty that needs to be dealt with.
To demonstrate how the control loop of the coupled simulation works, Figure 15 shows important simulation results in the power and gas grid for ten exemplary hours with a P2G device of 10 MW installed capacity. We compare three different cases: The base case which does not include any optimization or RES curtailment, thus leading to constraint violations in the power grid. The power grid restricted (PR) case includes RES curtailment and a P2G operation with the help of an OPF, but no feedback from the gas grid. Constraint violations in the gas grid might occur, as the P2G operation only follows the requirements of the power grid. The power and gas grid restricted (PGR) case also considers the gas grid constraints and the P2G power is reduced in case of violations, leading to a higher curtailment of RES plants. The left panels of Figure 15 show results of the power grid. The cumulated RES feed-in has to be reduced drastically in most time steps in order to comply with the power grid constraints. Otherwise line loading and voltage exceed the limits by far. On the right results of the gas grid are displayed. Here the feed-in from the P2G device has to be reduced in some time steps in order to comply with the gas grid constraints. The mass flow at the slack that represents the feed-in from the high pressure gas grid, would in some time steps become negative which does not comply with the uni-directional flow constraint. In some time steps also the maximum pressure at the connected gas grid node would be exceeded. Only in the PGR case all constraints are satisfied at all times which cannot be guaranteed in the PR or base case.

Study of the P2G Flexibility Potential
As the presented approach ensures a secure operation of both power and gas grid in the PGR case, it can be used to study how P2G devices can effectively be integrated into a coupled infrastructure. We performed a large scale analysis for different RES scenarios and different numbers, positions and sizes of P2G devices, leading to a total of 462 different configurations. The possible positions of the P2G devices with their connections to power and gas grid are shown in Figure 16. The option of direct hydrogen feed-in is only considered for position 1. In all other positions the fluid composition would not be the same everywhere, which is not possible to model with the applied version of pandapipes. An overview of the configuration setup is given in Table 7. Performing such a large study is only possible with the help of a parallelized approach conducted on a computer cluster, as an OPF is very time consuming. Relevant questions that can be addressed with the analysis are to what degree the P2G device can reduce RES curtailment and what influence the gas grid constraints have on this result for different configurations. As the main economic factor in this setting is the curtailed energy, we look at the results of a full year simulation. For the three RES scenarios Figure 17 shows results of curtailed energy from RES and the reduction of this curtailment through the P2G devices in relation to the installed P2G capacity (Relative Curtailment Reduction -RCR). The latter is an important key performance indicator, as installation costs of P2G devices mainly scale with the installed capacity. We look at four configurations, one without P2G device, one with the highest RES curtailment reduction, one with the highest impact of the gas grid constraints and one with the highest RCR. As can be seen, the potential to reduce curtailment is especially high in the case of the mix scenario. This is mainly due to the more distributed and lower RES feed-in peaks compared to the PV or wind scenario, which can more easily be absorbed by the P2G devices. The influence of the gas grid constraints is the highest in the PV scenario which is to be expected, as the highest PV peaks occur during summer, when gas consumption is usually very low and thus the gas feed-in is limited. Yet the impact is rather small: the P2G potential is reduced by less than 20% in all scenarios for the configurations with the highest curtailment reduction. In the configurations with the highest RCR this is even less, which raises the question whether a detailed gas grid simulation is necessary in this kind of analysis. It might be possible to identify simple rules that can be considered as constraints in the OPF. Furthermore it is obvious that position 2 is very promising for placing a P2G device, which is due to its sensitivity to line overloadings in the power grid near that position.
Another result of interest is how the potential to reduce curtailment changes with the installed capacity and number of P2G devices. Figure 18 shows the RCR for one, two and six installed P2G devices with capacities ranging from 1 MW to 10 MW in the PGR case. For different positions and RES scenarios, the results vary greatly, as can be seen from the box plots. In general, the curtailment is not reduced by much with increasing P2G capacity, so that the RCR decreases. There are mainly two reasons for this. On the one hand the higher P2G power is required less often to resolve all congestions    Figure 17. Analysis of the curtailed RES energy and its relative reduction by the P2G device in relation to the size (right-aligned due to the inverse behavior compared to the curtailment) for the three RES scenarios and different configurations. The simulated cases PR and PGR are shown in different colors to highlight the influence of the gas grid constraints on the results. Considering gas grid constraints leads to a lower P2G employment, so that more energy from RES needs to be curtailed. due to a lower occurrence of high RES peaks. On the other hand, the higher P2G power is more often restricted by the gas grid constraints. It also means that the RCR is expected to increase with further decreasing P2G size.
The previous analysis is of qualitative nature and it is not possible to derive the optimal P2G capacity from it. There are many unconsidered factors that influence whether a P2G device leading to a specific RCR can be operated economically. The cost of installing and operating a P2G device need to be compared to RES curtailment costs and grid expansion costs. Furthermore there might be other business cases for the P2G device. Thus it is necessary to compare different options, such as operating the device at times of low electricity prices, implementing a local flexibility market or direct contracting for delivering ancillary services. These options could even complement each other. Analyzing these questions would require an extension of the methods we presented.

Summary and Conclusion
We presented the new open source simulation tool pandapipes. When coupled with pandapower through the multi-energy module, it forms a powerful environment for the coupled simulation of MEGs. To the best of our knowledge, it is the first such simulation environment that is available open source and does not require proprietary software to work. After analyzing the field of MEG simulation approaches and tools, we defined three criteria that should be addressed by a dedicated tool: clear data structure, adaptable MEG model setup and performance. In the following we showed how pandapipes addresses these criteria in its implementation. Then we demonstrated two use cases: In the first we analyzed local peak shaving strategies in combination with district heating grids, in the second we performed a large study on the potential of a P2G device to serve as flexibility in a coupled power and gas grid. Although these use cases are based on artificial data, they are examples of typical research studies and show that the pandapipes tool is very powerful at performing coupled simulations and they underline the importance of the criteria with respect to large research studies.
We illustrated how the pandapipes architecture addresses the criterion of a clear data structure: • The use of pandas tables as input and result storage makes the structure very transparent and users can easily get an overview. Input data for models can easily be adapted manually or in an automated approach. Data from time series can be transferred to a new grid state performantly and for the result evaluation efficient and convenient pandas functions can be used. The similar design of pandapipes and pandapower facilitates getting started with the respective other tool and getting an overview over complex coupled infrastructures. • Also the internal structure of pandapipes is very clear. Component models are defined such that new models can be added without breaking the internal logic of the program. The use of numpy arrays with a unified design makes the internal structure clear for developers and parts of the calculation logic can be adapted without putting much effort into its restructuring. • The use of pandas also supports a simple setup for time series simulations, as the time series values can simply be connected to the input tables of the pandapipes grid via indexing. The same is true for the output which is related to the usual pandapipes output tables. • In the presented use cases, we used the pandas functionalities for post-processing. Even large amounts of data from time series simulations can be evaluated easily with respect to different characteristics that are of interest in the current research problem. For instance, we analyzed the full year time series results of more than 400 different simulation setups to identify how the size of a P2G device influences its potential to reduce RES curtailment. This analysis can be performed within seconds once the results are stored.
We showed how a multi-energy grid with different coupling models can easily be created with the help of the pandapipes multi-energy module, thus fulfilling the criterion of an adaptable MEG model setup: • The coupling of different infrastructures is done with the help of a multi-energy grid structure using multi-energy controllers. Models for coupling components can be chosen from the existing implementations or defined by the user himself. As controllers are implemented in an object-oriented approach, new controllers can be built upon base classes that define the simulation interactions. The user can focus on the physical model implementation. The controllers can implement arbitrary models of grid coupling components and can therefore be viewed as an implementation of the concept of energy hubs. • Once the multi-energy grid structure is built up, the user can start a time-series calculation.
During the simulation, all controllers that are used to couple grids included in the structure are called automatically. Information to be exchanged between the grids are transferred during run time in every time step, which is crucial, as the control algorithms of many coupling facilities require information from all connected grids for convergence. • The permissive open source license enables the adaptation of any parts of the software. Being based on python, many external libraries can be used to enhance modeling and reduce the work load for developers who want to extend the simulation framework. Python also allows console-like scripting, making it possible for unexperienced users to quickly learn how to use the simulation environment with just a few specific commands. • We showed the flexibility of the models in the use cases. In one case a heat pump controller was used to resolve line overloadings in a power grid. In the other case a gas grid controller was coupled with an OPF to ensure constraint compliance in the power and gas grid. This integration of different models within the tool structure is necessary to address a wide range of research questions.
We also performed a comparison between pandapipes and STANET®with respect to simulation time in order to show how pandapipes fulfills the performance criterion: • In comparison with STANET®, pandapipes was shown to be more performant. We distinguished simulation time and total time, where the total time includes a simulation overhead of input and output processes. For both simulation and total time we found pandapipes to be faster in every simulation setup. Reasons for this might be the more precise calculation, the slower readout speed as well as the use of a GUI in case of STANET®.
• Especially with respect to grid size, the simulation time does not scale as much in pandapipes as it does in STANET. It is therefore especially well suited for the simulation of complex grids that cannot be simplified with respect to the number of nodes and customers connected to them. • In [17] pandapower was also shown to be a very performant tool. As a coupled simulation can be performed within python and all data will be available throughout the whole simulation environment during the whole simulation process, this hardly induces any coupling overhead with respect to computation time. This is a large difference to a co-simulation approach, which requires a large overhead for data exchange and sometimes also for triggering single simulations externally. Co-simulations are rather well suited if further tools, such as market models, are required, or if exchanged data shall be limited due to confidentiality. • The performant core makes pandapipes a perfect tool for extensive investigations including e.g. probabilistic grid planning, time series, Monte Carlo and placement studies. This aspect is enhanced by the permissive open source license which supports running several simulations in parallel on several cores, a necessity for for many research studies. The demonstrated use cases underline how crucial performance is, as a large number of simulations had to be performed for each of the studies.
Although the presented simulation environment covers a large variety of use cases, there are still model limitations that shall be addressed in future releases of pandapipes. For instance, up to now the heat transfer option is restricted to incompressible fluids. Although a user can do a heat transfer calculation using a compressible fluid, the results may not be accurate enough, because in pandapipes the calculated temperature does not influence the hydraulic fluid properties like density or viscosity. Transient calculations can be decisive in modeling effects in district heating grids, especially when integrating heat storages. Another missing feature is the calculation of fluid mixtures and their propagation in a piping grid, especially for the evaluation of use cases where large amounts of hydrogen are introduced into a natural gas grid.
Other extensions shall increase the range of functions and the performance in order to make pandapipes more attractive for specific use cases. They include a speed-up of time series simulations, new or extended libraries for component models, controllers, topology analyses and grids for testing purpose. With the increasing interdependency of infrastructures also an integration of power and piping grid simulations on kernel level will become more important. This is especially true for OPF calculations, where an optimal operation strategy must comply with the constraints of all integrated grids. With such an integrated OPF, the flexibility provision model could be extended so that the optimization is performed in one step without feedback loop. For this purpose an approach similar to [12] or [65] could be chosen, or a port to existing tools could be used, like the GasPowerModels Julia package [73].

Appendix A. Equations describing hydraulic properties
For an incompressible medium, the pressure variation along a pipe section with length dl can be described using equation A1. The pressure difference is caused by three effects, described with the terms on the right-hand side of the equation: The pressure loss due to friction is proportional to the darcy friction factor λ. It also depends on the pipe diameter d.
A lumped pressure loss can be expressed using the pressure loss coefficient ζ. This pressure loss is introduced e.g. by bends or components which are not represented by a component model. Both the pressure losses due to friction and to lumped components are proportional to the dynamic pressure ρ·v 2 2 , where ρ is the fluid density and v is the flow velocity.
A pressure difference introduced by a given height difference ∆h between the connected junctions is given by the last term. This difference depends on the fluid density and the gravitational acceleration g.
In case of compressible media, the density of the medium is expressed by using the ideal gas law and a reference state, as shown in equation A3.
The reference state, noted with the index N, corresponds to the known gas normal properties.To add real gas behaviour, the equations resulting from the description as an ideal gas can be modified by a compressibility factor.
If relation A4 is used to express the relation between the mass flows of the system and reference states and the density ρ is additionally replaced by expression A3, the system velocity can be formulated as a function of the velocity of the reference state v N (equation A5): Inserting equations (A4) and (A3) into (A1), leads to equation (A6) (including only pressure losses due to pipe friction): This makes it possible to determine the pressure losses in the pipe based on the completely known reference state. Because this equation only contains the gas velocity with respect to the reference state c N , the gas velocity for the system has to be determined after the calculation of pressure losses from equation A4.
Compressible media lead to a continuous velocity and pressure distribution along a pipe. To capture these distributions, derivatives with respect to the pipe coordinate of equation A6 are described with the finite difference method. As a result, the user can divide a pipe into several subsections, gaining a finer resolution of solution variables for the pipe of interest.

Appendix B. Equations describing thermal pipe properties
If the heat transfer mode is started, equation A7 is solved for a pipe section with length dl. To solve this equation, the mass flowṁ has to be known. It can be determined by a preceding hydraulic calculation. Because equation A7 does not contain a partial derivative with respect to time, temperature distributions do not include inertial effects.
In case of equation A7, the right-hand side describes a heat flow exchanged with the ambient, which has a temperature of T ext . This exchange is described via the heat transfer coefficient α. Besides, a user-defined optional heat flow can be defined using the term q ext .

Appendix C. Model Implementation in pandapipes multi-energy
When modeling coupling facilities, a multitude of different control schemes can be chosen. In Figure A1 controllers for a heat pump and for a P2G are illustrated with the information that is exchanged between the controllers and the different grids. This data exchange includes the input and output of control and reference variables and at the same time ensures a correct model of the energy transfer. The control loop with checking controller convergence in combinations with pipeflow and powerflow calculations is built around these models.

Appendix E. Additional information on the peak shaving scenario
This section provides some additional information about the heat pump and prosumer controller. The control strategies implemented for the heater components and the heat pump are shown on Figure A5.
Both algorithms introduce two set points, a primary (T_set_pr) and a secondary (T_set_sec), to control temperature values inside the connected storages. The left side illustrates the control for the electrical heaters. If the storage temperature lies above the secondary setpoint, the heater is switched off, because the maximum temperature of the storage is exceeded. Otherwise, the algorithm checks if there is enough excess power provided by the PV plant to switch the heater on. If this is not the case, the heater is still switched on if the temperature is below the primary set point minus an offset of 2.5 K. In this case, the necessary power is provided by the power grid. Assuming that there is still not enough excess power provided, the heater is switched off if the storage temperature is above the primary set point plus an offset of 2.5 K. If the temperature is between the offsets, the current state of the heater remains unchanged. This introduces a hysteresis in order to avoid a rapid change of state of the heater.
The control algorithm of the heat pump works in a comparable way. The most important difference is that the algorithm checks for a line overloading if excess power is available.   Figure A5. Implemented control algorithms for the electrical heaters (left) and the heat pump (right).

Acknowledgments:
The authors thank Christian Spalthoff (Fraunhofer IEE) for his design and illustration support.

Conflicts of Interest:
The authors declare no conflict of interest.

Abbreviations
The following symbols and abbreviations are used in this manuscript: Preprints (www.preprints.org) | NOT PEER-REVIEWED | Posted: 6 October 2020 doi:10.20944/preprints202010.0132.v1 Table 6. Constraints and boundary conditions for the grid simulations. Voltage and current limits apply for the whole power grid. Pressure and hydrogen fraction limits apply only for the power-to-gas (P2G) connection node, as the respective maximum value will always occur there. N HH stands for the number of households connected to the gas grid.
power grid gas grid preset voltage at slack node V slack = 1.01 p.u. preset pressure at slack node p slack = 0.7 bar maximum voltage V max = 1.06 p.u. maximum pressure P2G node p P2G,max = 0.75 bar maximum line loading I max = 0.6 I max,th maximum hydrogen fraction x H 2 ,vol,max = 10% uni-directional flow at slackṁ max,P2G = ∑ N HH l=0ṁ l