Modeling multimodal energy systems

: Information and communication technology (ICT) and the technology of coupling points including power-to-gas (PtG), power-to-heat (PtH) and combined heat and power (CHP) reshape future energy systems fundamentally. To study the resulting multimodal smart energy system, a proposed method is to separate the behavior of the component layer from the control layer. The component layer includes pipelines, power-lines, generators, loads, coupling points and generally all components through which energy flows. In the work at hand, a model is presented to analyze the operational behavior of the component layer. The modeling problem is formulated as state and phase transition functions, which present the ex-ternal commands and internal dynamics of system. Phase transition functions are approximated by ordinary differential equations, which are solved with integral methods. State transition functions are nonlinear algebraic functions, which are solved numerically and iteratively with a modified Newton–Raphson method. In a proof-of-concept case study, a scenario shows the expected multi-sector effects based on evaluated models.


Introduction
gradually by experts from different disciplines with different goals, requirements and limitations. In addition, coupling points which interlink different energy sectors make the resulting CPES more heterogeneous and deviant from the original specification. The controllers of these systems could not be designed and parametrized uniformly. The reason for this lies in the nonlinear and time variant decision making process of autonomous controllers. An illustrative example is given in [24] where different reaction times of transformer controllers and solar cell inverter controllers lead to a voltage fluctuation.
A specific class of controller interactions that can be understood as unwanted emergence of a complex system, are called controller conflicts. In [18], a framework with three layers to identify these conflicts is presented (see Fig. 1). First layer of the framework is modelling MES, which involves the integrated models of the energy sectors and coupling points. This model has to satisfy the two following requirements: 1. The model of the MES must have a high degree of abstraction. The resulting model should show the mutual interaction between different energy sectors and includes coupling points. The model should reflect the states of energy systems by changing energy balance. 2. The related control system of the MES should not be embedded in the model to allow for a later exchange and evaluation of different control concepts. Figure 1: Conceptual view on the targeted framework for controller conflict analysis in multimodal energy systems [18]. In the work at hand, the formalization of the lowest layer, i. e., the layer of physical component simulation, is presented.
The resulting model serves as an environment in which controller actions, interactions and conflicts can be presented and evaluated. The demanded degree of abstraction should be chosen as the functional behavior of MES. Within this study, the multimodal concept is restricted to power, gas and heat at distribution level as shown in Fig. 2.
The coupling points -including methanation (PtG), heat pumps (PtH) and gas engines CHP -are regarded as active players in controlling the energy system. The MES model is assumed to be enveloped in an energy balance zone, although gas and electricity networks have connections to the transmission level. The gas sector is composed of networks of pipelines and components including compressors, valves and gas regulators responsible for transportation, distribution and delivering gas with a specific pressure and quality to endusers [1]. In national gas systems, the gas is compressed through compressor stations to a defined pressure level and is kept at this level. Through gate stations, gas flows from national gas distribution networks to the regional gas distribution networks, which are modeled in this study. These networks do not include compressors and show a pressure range of 1 to 7 bar.
District heating networks are regional networks of pipelines and components including circulating pumps, pressure control systems and heat exchangers, which deliver heat to end-users within certain pressure and temperature ranges. These networks are typically connected to a power plant as the main heat resource; in case, the supplied heat is exceeded by heat demand, the remaining demand is compensated by other heat generators such as heat pumps as considered in this study.
Electrical power grids include power transmission and power distribution systems. The focus of this research is on distribution systems at high and medium voltage levels. In addition to passive networks of power-lines, inverter-based generators, synchronous generators and local power transformers are modeled in this study.
This contribution presents the model development process to an executable MES model with considering the coupling points as active controlling components. Thus, the presented model serves as the basis for the development of adaptive control concepts in later work. In section 2, the state of the art regarding MES modeling is presented to illustrate the knowledge gap in modeling MES. In section 3, the modeling, its preliminaries, the related assumptions and the implementation of the model are explained. A case study is introduced in section 4 as proofof-concept and the results are discussed. In section 5, the work is concluded with an outlook on future work.

Related work
MES modeling is discussed using the terms sector coupling, integrated energy systems and multi-energy systems. The studies can be categorized into planning, operational planning and operation of MES. The planning of MES often includes finding the position, dimension and capacity of components, e. g., finding the dimension of PtH in comparison to hot water storage tank in a district heating system [20]. In [14], an overview on the concept of MES is presented, including the different models from a techno-economic perspective. In works related to operational planning [5,14], the role of coupling point technologies for improving the energy supply performance and enabling economic dispatch is discussed. The main goal of these works is to find the optimal scenario for energy flows between different energy sectors. One of the proposed solutions has been conceptualized as energy hub in [8]. None of these studies aims to analyze the MES controllers or the MES state. For instance, they do not consider the effect of heat demand on the gas pressure. All these works have a higher abstraction degree in comparison to the presented work. Our modeling focus is on operation of MES where the reaction of controllers and the behavior of MES is the main topic, viz., we consider more details in comparison to the above mentioned works to allow for a controller conflict analysis in later steps.
In [16], possible applications of coupling points in electrical power grids (transmission and distribution) are discussed. The study points out the provision of ancillary services by coupling points. The authors discuss frequency control and the voltage control as examples for ancillary services. Our study assumes that the coupling points have the potential to provide ancillary services bidirectionally. Therefore, we use this assumption for modeling the MES. Some other works have studied the side effects of sector coupling and discussed the operational behavior of every energy sector under consideration of coupling points. The effect of feeding gas by PtG into gas networks is analyzed in [1]. The results of this works are deployed to choose proper coupling point technologies. In [15,13], integrated gas and electricity networks are investigated. The proposed models consider merely the energy flow from gas sector into electricity sector. They have a lower level of abstraction and consider more details. Additionally, a separation of control and component level is not given.
In [12], the authors sketch out an integrated power grid and district heat with considering CHP as a single unidirectional coupling point. The authors of [12,21] have developed an integrated static model for gas, heat and electricity sectors, which covers the system behavior at discrete time points. We have developed this approach in relation to our requirements. Our model has to be able to point out the MES behavior in stationary state and between two consecutive stationary states as well. In addition, this study models distributed slacks for all networks, which facilitates modeling an energy system with different connections between sectors.

Methodology
The steps of model development in this study are as follows: 1. Determining the operation of MES 2. Determining the states of MES 3. Discussing the time response 4. Modeling external supply 5. Modeling external demand 6. Modeling networks 7. Modeling coupling points 8. Model development 9. Implementation and simulation.
The first step determines the operation of MES in form of energy balance equations. With respect to the operation of MES, the state of MES is defined in the second step. Since the state of MES originates from different energy sectors with different time constants their temporal behavior is discussed in the third step. Steps 4 up to 7 involve modeling the components of MES viz., finding mathematical equations which relate energy balance equations to the state of MES. In step 8, all component models are used to construct an integrated mathematical model for MES, and step 9 comprises details regarding the implementation and simulation of the mathematical model. 1

Functionality of the model
Energy supply systems have to ensure the fulfillment of energy demand with a specific quality, i. e., the main duty of energy system operators is to keep the energy system stable and within the defined operational limits. Stability of an energy system across all sectors means e. g., that the deviation between supply and demand must be always fewer than a specified value. The energy demand in course of time (i. e., load profile) is one input to the model. The energy supply should be controlled in a way that supply follows the demand. Hence, the supply energy, which is determined by the control system, is another input to the model. 2 Obviously, there is almost always a mismatch between these two inputs. The proposed model will demonstrate the impact of these mismatches on certain nodes of the energy systems. The nodes of MES are normally physical points with specific characteristics, for instance, they are heat exchangers, gas connections, electrical bus bars or even a specific point at district heat whose temperature is of interest. In this work, we use the superscripts e, g and h to label electricity, gas and heat respectively.
The energy imbalance at nodes is modeled based on the conservation of mass and energy, and Kirchhoff's rules. If node i is part of the electrical power grid, the algebraic sum of inflows and outflows of active power p and reactive power q at this node should be equal to zero, as formulated in equations (1) and (2).
1 As in this contribution a proof-of-concept is presented, no detailed scenario setup is presented for reasons of brevity. For future studies, a refined approach will be used following the design of experiments approach.
2 Please note that this is a modeling assumption interacting with ICT-based control that will be analyzed in later steps. While supply and demand can be dynamically adapted using various ICT-based approaches, in the modeling concept, the balance has to be realized.
In case that node i is part of the gas sector, the conservation of mass implies that the sum of inflows and outflows of mass flow rateṁ at point i is equal to zero. Considering a constant density of gas, mass flow rate is converted to volumetric flow ratev . Equation (3) points out the balance for gas system nodes.
To maintain conformity of the model, gas demand p g with unit J s is converted using equation (4) into volumetric flow rate where GCV stands for gross calorific value of gas [1].
If the node i is part of the district heat sector, the conservation of mass and energy implies equations (5) and (6), in whichṁ stands for mass flow rate and u represents thermal energy flow.

Determination of states
In the work at hand, the state of energy system is represented by control variables which specify the quality of delivered energy to consumers, and energy system operators are bounded by contract to keep them in a predefined tolerance band. Technically seen, energy imbalances can only be identified by measuring these variables. For example, not any active power deficit in the electricity sector can be measured directly, but can be observed indirectly by measuring the network frequency. In the electricity sector, the amplitudes and phases of node voltages v and θ and the network frequency f as well are control variables. The node pressure ψ g is the control variable of the gas network. Relevant control variables for the state of the district heat system are node pressure ψ h and node temperature ϑ.
The variables representing the state of MES nodes are positioned in the vector X and set up the state vector of MES as shown in equation (8). The goal of modeling MES is to find the value of vector X at every arbitrary time t.

Time response
At every point in time the MES is either in a stationary state or in transient state. In stationary state, vector X is timeinvariant, which is addressed with X j . In transient state, vector X is time-variant, shown by X(t). Starting from a stationary state X j , any change in energy balance leads to a transient behavior X j (t). If we assume that an energy system has BIBO 4 -stability, this transient state ends up in the next stationary state X j+1 . We can deal with the transient behavior with different aspects. Many agent based studies assume that the time response of the energy system is fast enough and that it reaches its next stationary state timely, therefore there is no need to pay special attention to transient behavior. As a result, the system can be assumed to be always in stationary state and algebraic equations represent the energy system behavior properly. This assumption is not correct in the case of MES: firstly, the change rate of demand is high and in many cases higher than the energy system response. Secondly, the MES includes different sectors with different reaction times and different time constants, which have a range from milliseconds in the electricity sector up to hours in the gas sector. As a result, the stationary as well as the transient behavior of system state has to be considered. The former will be modeled with state transition and the latter with phase transition functions. As depicted in Fig. 3, the model consists of algebraic equations, which represent the state transition of the system and differential equations, which pinpoint the phase transition of system. State transition functions describe the system statically by considering external stimulation. The phase transition functions consider the internal dynamics of the system, which are modeled as first and second order transfer functions as suggested in [11] and [3] for the gas sector.
In what follows, energy system components are described using four categories: (a) external supply, (b) external demand, (c) networks, and (d) coupling points. The

External supply
External supplies in MES include import of gas from national gas networks to the local gas networks, feed-in electricity from overlaid networks, e. g., from high voltage power grid to the medium voltage power grid, and heat supply from power plants. These external supplies are connected to the so called slack nodes. As an example, the regional gas system with a pressure range from 0.7 to 7 bar is connected via a gate station to the national gas system with the pressure of 70 bar. The gate station is responsible to keep the pressure at regional side constant. Slack nodes are assumed to be reference point with known values in vector X, which are infinite sources of energy. Thus, slack nodes have known complex voltage value in the electricity sector, known pressure in gas sector and known pressure and temperature in district heat as shown in Tab. 1.

External demand
External demands in MES are gas, heat and power consumers, which are connected to load nodes in the specific sectors. In contrast to slack nodes, their energy demand and load profile is assumed to be known. Therefore, load nodes have fixed values in the vector F, whereas their states are unknown. The external demands are modeled as passive loads which means that the drawn energy does not depend on the energy system state. Strictly speaking, load nodes have known active and reactive power, gas demand and heat demand depending on their sectors as presented in Tab. 1.

Networks
MES include networks of power lines, gas pipelines and heat pipelines for transporting energy and distributing it to consumers. These networks include at least two groups of nodes: slack nodes and load nodes. In addition, any arbitrary point of MES which is necessary to be observed or controlled can be added to the node set e. g., junctions (intersections) in the gas sector and both sides of transformers in the electricity sector. These networks have been formulated according to the definition and methods from graph as described in [17]. Any graph is defined basically with its node and branch sets. 5 Based on the definition of nodes, the branches of MES are defined as any connection between two nodes with energy flows. The MES branches are discussed later in detail. Based on the defined nodes and branches of the MES, the incidence matrix is generated, which facilitates modeling of the network. In practice, the incidence matrix paves the way to relate energy balance vector F to the vector X by means of governing equations of the branches.
In electricity sector, a branch can be a power line or a transformer between nodes i and l. They are modeled in this work as 4-pole π-models with a series admittance y s and two identical parallel admittance y p . The resulting admittance matrix y il is shown in equation (9). The incidence matrix must be adapted for 4-pole models as well [6].
This admittance matrix is used to calculate the complex value apparent power s il = p il + jq il which flows between nodes i and l. It is calculated with equation (10) and is part of vector function F. In this equation, v i and v l represent the complex value voltages of nodes i and l, which are part of state vector X. 6 The result of the presented modeling approach are nonlinear algebraic equations for calculating the vector X at time point t. It has been assumed that the network of power lines and transformers has relatively small storage capacity and consequently very small dynamics, which can be neglected. 7 Thus the electricity network cannot change energy system states between time points t 1 and t 2 on its own.
In district heat and gas sectors, the branch can be a pipeline which connects node i and l. To avoid repetition, only gas equation will be addressed in detail in the following. The temperature in district heat network will be modeled separately. The pressure drop between nodes is owing to the flow resistance in pipelines and is modeled in this work with Polyflo equation, which is a simplified version of the general gas flow equation (Weymouth equation). The Polyflo equation prepares a proper approximation for the networks with a pressure range within [0.75, 7] bar. The equation (11) shows the pressure drop between two nodes i and l where L, D and E stand for the length, diameter and efficiency factor respectively. The pressures ψ i and ψ l are members of vector X and the gas volumetric flow ratev il is used in the vector F [19].
The same process is carried out for in-compressible flow in pipelines of district heat sector according to the Darcy-Weisbach equation. The results are nonlinear algebraic equations, which relate the node pressure of heat sector ψ h , which are elements of vector X, to the mass flow rateṁ, which are parts of the vector F. The internal dynamics of a pipeline between two nodes i and l and between two time points t 1 and t 2 is modeled with the time varying pressure function ψ g (t). Partial differential equations are used for analyzing the dynamics of pipelines as shown in [19]. The authors have suggested in [10] and [3] that the dynamic behavior of pipelines can be estimated as second order ordinary differential equations with an acceptable deviation. Equation (12) shows this relation between incomingv in and outgoing ψ g out with Laplace operator s and the transfer function H g . The constants in equation (12) are dependent on pipeline parameters and taken from [3].
The node temperature ϑ in the heat sector is another element of state vector X. For thermal calculation, every pipeline is assumed as a heat exchanger which dissipates thermal energy into its environment. The thermal model of the pipeline between nodes i and l in stationary state is described in equations (13) and (14). These equations relate elements from vectors X and F.
where z is the correction factor and is assumed to be equal to 1, χ is the radial heat transmission coefficient, 8 πdl is the surface area of pipeline and Δϑ lm is the logarithmic middle value of the temperature between environment and the working fluid [23]. Besides the states including pressure, temperature and voltage, which are defined subject to nodes, network frequency f is a global state variable, which is assumed to have only one value for all the whole electricity sector. The frequency of a network with only one synchronous generator is determined via the rotational speed of synchronous generators as shown in equation (15).
Parameter m is the inertia moment of the generator, p mech is the mechanical input power into the generator, p elec is the outgoing electrical power which is fed in the network and ω refers to angular frequency. 9 Parameter τ is a positive real value representing the droop control of the generator. Furthermore, τ considers the internal friction of the machine, which acts as damper [7,25]. If the power grid has g generators including synchronous generators and inverter-based generators like solar systems, and if the network is in stationary state, the left side of equation (15) for all generators is equal to zero and the network frequency is equal to f net . Substituting f net in all equations and summing them up, the network frequency is calculated by equation (16).

Coupling points
Three types of coupling points, PtG, PtH and gas engine CHP are discussed in this part. Coupling points can generally be modeled in two ways: (1) If they are modeled as two nodes, they are a load in one sector and generation in the other sector. (2) Coupling points can be modeled as branches as well, which relate energy balance of one energy sector to energy system states of other sectors. The latter approach considers coupling points as components which shift power imbalance from one sector to another sector. This approach is deployed in this study. As a result, coupling points have a direct impact on power balance and an indirect effect on energy system states. By analogy, coupling points as branches have the same function as gyrators in the bond graph modeling method. The coupling point branch is between two nodes. The inflow side is defined as load node, and the outflow side as slack node. The details of the modelled coupling points including functional behavior and boundary conditions are presented in the upcoming sections.

Power-to-Gas
PtG technologies include an electrolysis process which generates hydrogen and a methanation process which produces methane. The allowed amount of hydrogen in a gas network is limited to about 1 % of the capacity of the network [1]. In addition, hydrogen has a significant lower gross caloric value in comparison to methane. 10 Hence, considering hydrogen as output of PtG unit needs taking extra actions. If we define that PtG units produce only methane, then the gas mixture in a gas network is homogeneous. In a PtG unit, water, carbon dioxide and electricity are the inputs and methane, oxygen and water are the outputs of the process. The backbone of the PtG process consists of electrolysis and methanation in a row as shown in Fig. 4. For both processes, electrolysis and methanation, there are some different technologies. The modeling of PtG in this study needs a high degree of flexibility including short reaction time and a wide operation domain from offstate up to rated working point. At present, these requirements are nearly met with alkaline electrolysis (AEL) technology. A load changing from 5 % to 100 % and a response time of 0.5 seconds is possible [9]. Without limiting the generality, this technology is chosen in course of this work. For methanation process, we assume that there is a source for carbon dioxide with the needed purity of carbon dioxide. The methanation process can be done with two different technological ways: the biological and chemical methanation. The former is a slow and inertial biological process, the latter is in comparison a fast chemical thermodynamic process. The chemical process converts hydrogen and carbon dioxide as low density energy carriers into methane as a high density energy carrier. This exothermic reaction has an efficiency of about 83 % and the product is delivered with a low temperature and high pressure. Thus, the chemical process is chosen and modeled in this work. For calculating the time constant of PtG unit, the plant volume and final gas velocity are important parameters. The PtG unit which is modeled in this work has a minimum load of 20 % of rated load and the efficiency is assumed to be 50 % [9].
The PtG units are modeled as a branch which has a slack node in gas sector and a load node in electricity sector. The slack node injects the mass flow rateṁ in the gas sector and demands active power p from electricity sector as shown in equation (17), where η PtG and GCV stand for overall efficiency of PtG process and gross calorific value respectively.ṁ

Power-to-Heat
Without loss of generality, electrical heat pump (HP) is considered as a chosen technology for PtH units. A decentralized HP is modeled with a constant source temperature (T source ), controllable pressure and constant coefficient of performance (COP). The minimum load is considered to be 50 % of the rated load. The modeled HP operates at following operation steps: 50 %, 60 %, 70 %, 80 %, 90 % and 100 % of rated heating capacity. The dynamic behavior of the heat pump is modeled with a first order differential equation. Similar to the case of PtG, the PtH units are modeled as branches, which have a load node in electrical sector and a slack node in heat sector. The relation between active power p and the generated heat u is sketched out in equation (18). The HP is modeled to be a slack node in heat sector, therefore the generated heat is fed in with constant pressure and constant temperature [23].

Gas engine CHP
The gas engine CHP units have impact on all three sectors at once. In contrast to PtG and PtH units, they generate electricity, therefore they shift energy imbalance in the opposite direction as PtG and PtH. The thermal power from gas combustion is converted first to mechanical power p mech and via a synchronous machine is converted to electrical power. This conversion impacts the network frequency according to equation (16) and the internal frequency of synchronous machines according to equation (15). The gas engine CHP units operate either strictly electricity-or heat-driven. The electricity-driven CHP units are modeled as two branches with a load node in gas sector, a slack node in electricity sector and a negative load in heat sector. The heat-driven CHP units can be modeled in a similar way. If CHP units would be tooled up with heat storage systems or electrical boilers and even combined with HP, then both nodes in electricity and heat sectors can be modeled as slack nodes as chosen in this work.

Model development
Up to now, the balance vector F, the energy system state vector X and nodes and branches of the energy system are defined. The behavior of coupling points are modeled as well. The entries if vector X and inputs and outputs of coupling points are depicted schematically in Fig. 5. The relation between balance vector F in all nodes of MES with state vector X has been described based on the modeling of individual components of the network. The resulting vector F(X) includes nonlinear algebraic equations which determine the behavior of energy system states in stationary state. To solve these equations normally no analytic solution can be applied, but rather they are solved numerically and iteratively. In this study, we have employed Newton-Raphson method and equip it with an extra calculation for estimating initial values to improve the convergence of calculation. To determine the behavior of en-ergy system states between every two time points, the ordinary differential equations of components are solved with the integral method (Fourth-order Runge-Kutta) [7].

Implementation and simulation
The proposed MES model is implemented in MATLAB. A general class for MES is defined. Properties of this class include network structures and branch characteristics. Properties of this class read data from Excel-sheets, convert values in per unit, build the graphs and their incidence matrices, finding suitable initial values for vector X, perform load flow calculations based on Newton-Raphson method, calculate the values of vector X in stationary and transient states and visualize the results. Every object of this class represents a specific MES. The flow chart of load flow calculation is discussed in earlier studies [1], but the methods to find initial values for vector X and running load flow calculation for MES are the most challenging steps in the implementation. The load flow equation of gas and heat (equations (5) and (3)) are s-shaped, which make Newton-Raphson method sensible to initial values. As a result, the initial values have to be found with special care. Concatenating states from different sectors in one vector leads to convergence problems in Newton-Raphson method as well. To overcome this problem, the Newton-Raphson method is tooled up with Bisection method [4].

Case study
The presented modeling for MES is tested as a proof-ofconcept case study to show the general feasibility of the approach. The case study includes networks from all three sectors. For the electricity sector, a meshed test system with 9 nodes is chosen. As an exemplary gas system, a model with 11 nodes from [1] is used. The test system of heat network is a simple district heat network with 9 nodes and 8 branches, which has a line topology. The test systems are interlinked with the PtG, PtH and gas engine CHP units as shown schematically in Fig. 6. The sectors heat, gas and electricity are labeled with colors red, yellow and green and circle-, diamond-and square-shaped respectively. The oversized nodes represent the coupling points between sectors.
For the case study, the energy balance vector F has 45 equations including 9 equations for active power, 9 equations for reactive power, 11 equations for gas volume flow, Figure 6: The graph of a proof-of-concept case study for MES modeling with three coupling points including PtG, PtH and gas engine CHP. The oversized nodes represent nodes of coupling points. 8 equations for mass flow of heat network and 8 equations for determining temperature of nodes. It is assumed that the gas sector has only one slack point at node 1, the heat sector has two slack points at nodes 1 and 9, and the electricity sector has 3 slack points at nodes 1, 2 and 3. The gas engine CHP, which is connected to node 1 of the electricity grid is a load at node 11 of the gas network and feeds into node 9 of the heat network. Nodes 6 and 8 of the electricity sector feed in the gas and heat sector via PtG and PtH units. The MES model is validated by comparing the results with the individual case studies by means of disconnecting networks and testing load flow calculation of single sectors. In the case of electricity, the network is configured with the same parameters from literature [2] and results have shown the validity of model. The results from gas system are compared with the results from [1]. The heat sector has the line topology and the results are calculated and compared manually. To show the functionality of the model and to test the fundamental requirements as depicted in section 1, a test scenario has been designed. In this scenario, the electrical load is kept constant and the gas demand is dropped off in course of time. The results are shown in Fig. 7. The figure shows the node pressures of the gas sector and voltage amplitudes of the electricity sector in course of time.
As expected, the gas pressure grows owing to smaller mass flow rate. Decreasing gas demand has the side effect on the electrical network, where the voltage level rises. All of the executed tests have proven that the goal of modeling is achieved, i. e., deficit or surplus of energy in one sector affects other sectors inherently.

Conclusion
The goal of the work presented here was to set up a formal model and executable simulation model for MES. This model acts as a component simulation layer for the analysis of controller conflicts in MES as depicted in Fig. 1. Based on the fundamental laws for conservation of energy and matter at every node of MES, balance equations have been defined. Furthermore, the state of MES is presented by operating variables of MES. The nodes are categorized as external demand, external supply and junction nodes. The network of nodes -representing the coupled energy systems -is presented as graph and matrices, which facilitate the implementation and make the modeling scalable. The time response of MES is divided into the stationary state and the transient state. With all these definitions, the next step is to determine the behavior of individual components based on the elements of state vector and balance vector. The resulting model has been implemented and tested in a proof-of-concept case study.
The proposed model includes some assumptions, which have to be studied in more details. For instance, the impact of lower GCV of injected synthesized gases to gas networks can be added to the model. Coupling points have been modeled as components which can shift energy from one sector to another. The implementation techniques have to be enhanced for a guaranteed convergence of Newton-Raphson method as well. In case of networks with a mesh topology, the initial estimation for Newton-Raphson method determines whether the MES states can be calcu-lated or not. The differential equations which represent the transient behavior of MES states are assumed as first and second order, therefore they can only sketch out a group of dynamical behavior. The related time constants have to be adjusted in relation to the sampling time of the control system that will be implemented as an agent-based system in the next step. We have followed the approach to model MES separated from the control system, therefore the next step of developing MES model is to determine the role of control systems to identify and mitigate controller conflicts.