Power-Transfer-Distribution-Factor-Based Sensitivity Factors for Integrated Energy Systems

The stronger the electric power, district heating and gas system are linked in an integrated energy system (IES) the larger their flexibility potential but also their complexity. To ensure a secure system operation, the effect of a coupling unit's power change on the overall system state must be analyzed before the power change is carried out. Otherwise, coupling units can unintentionally cause limit violations in one energy system. Currently, for such analyses a power flow calculation must be performed, which is computational expensive and time-consuming. This article presents sensitivity factors extending Power Transfer Distribution Factors from electric power system analysis to IES. The sensitivity factors include the dynamic behavior of district heating and gas systems as we derive the sensitivity factors from a joined quasi-steady-state power flow calculation. In two case studies we profoundly test the applicability of the IES sensitivity factors, their accuracy, and limitations. We show that the sensitivity factors are on average ten times faster in estimating a new system state after a unit's power change compared to a power flow calculation. Also, the sensitivity factors can provide good estimates and have a higher accuracy than sensitivity factors which are derived from a steady-state power flow calculation.

Power-Transfer-Distribution-Factor-Based Sensitivity Factors for Integrated Energy Systems Jonte Dancker and Martin Wolter , Senior Member, IEEE Abstract-The stronger the electric power, district heating and gas system are linked in an integrated energy system (IES) the larger their flexibility potential but also their complexity.To ensure a secure system operation, the effect of a coupling unit's power change on the overall system state must be analyzed before the power change is carried out.Otherwise, coupling units can unintentionally cause limit violations in one energy system.Currently, for such analyses a power flow calculation must be performed, which is computational expensive and time-consuming.This article presents sensitivity factors extending Power Transfer Distribution Factors from electric power system analysis to IES.The sensitivity factors include the dynamic behavior of district heating and gas systems as we derive the sensitivity factors from a joined quasi-steadystate power flow calculation.In two case studies we profoundly test the applicability of the IES sensitivity factors, their accuracy, and limitations.We show that the sensitivity factors are on average ten times faster in estimating a new system state after a unit's power change compared to a power flow calculation.Also, the sensitivity factors can provide good estimates and have a higher accuracy than sensitivity factors which are derived from a steady-state power flow calculation.Velocity in m /s.

I. INTRODUCTION
I NTEGRATED energy systems (IESs) can provide flexibility for the electric power system (EPS) and support the decarbonization of other consumption sectors.They can improve the system operation [1] and economic efficiency of the overall energy system [2].For example, when electrical transmission lines are congested due to a high generation of volatile renewable energy sources (RES), the generated energy can be shifted instead of curtailing the generation units (e. g., [3], [4]).Moreover, district heating systems (DHSs) and gas systems (GSs) can provide a considerable flexibility due to the thermal storage in the DHS [1] and the gas compressibility and potential of hydrogen blending in the GS [5].This does not only reduce cost in the EPS but also increases the share of RES in other energy consumption sectors.
The potential of IESs is higher the tighter the different energy systems are connected, i. e., the more coupling units are used.However, the interactions lead to a higher complexity in the planning and operation of IES [1].An uncoordinated operation of the coupling units can result in unintentional limit violations and pose a threat to the IES' system reliability [6].Hence, the effect of a power change of a single or multiple coupling units on the overall IES' system state must be analyzed before the power changes are carried out, identifying if a coupling unit can shift energy between energy systems without causing a limit violation.To ensure a secure IES operation while considering the whole flexibility potential, such analysis must include the system behavior and the interactions between the different energy systems [7].Although this analysis becomes more important and complex the more coupling units are installed, currently only power flow calculation methods allow such analysis.These however, are computationally expensive and are not suitable if many scenarios must be investigated, e. g., for congestion management or grid planning purposes.
For such analyses, sensitivity factors can be used as they describe the effect of power changes of multiple units on the power flow by linearizing the non-linear power equations [8].Once these sensitivity factors are derived, many different scenarios can be easily simulated without performing another power flow calculation, increasing the computational efficiency.Such sensitivity factors, including Power Transfer Distribution Factors (PTDFs), are widely used in EPSs as they have a high computational efficiency and accuracy [9].
Despite the great potential of sensitivity factors only few publications apply them to IESs.[10] and [11] use sensitivity factors for electricity-gas IESs.The approach of [10] is restricted to meshed networks and to the assumption that flows and pressures change instantaneously.Moreover, the sensitivity factors are determined iteratively in an optimal power flow optimization [10].The approach of [11], on the other hand, only determines the sensitivity factors for nodes at which gasfired generation units are placed and only determine the impact on the GS.[12] uses sensitivity factors in an optimal power flow optimization of a gas-heating-electricity IES.The sensitivity factors, however, only determine the effect of a wind power change on the operation of coupling units, aiming to adapt the operation of the coupling units to the wind power generation.
Overall, the sensitivity factors in the above-mentioned literature do not estimate the effect of a power change on the power flows in each energy system of an IES but only the effect on selected units.Thus, their factors do not include feedback effects which are caused by changing power flows.For example, a change in one system causes a change in another system which then affects the first system.Also, these methods only allow to investigate the change of one unit, which is not adequate in a tightly coupled IES in which many coupling units adapt their operation simultaneously.Moreover, the above-mentioned literature does not present a detailed analysis of the accuracy and limitations of the sensitivity factors, which is necessary to prove their benefits in possible use cases.In contrast, [13] presents PTDF-based sensitivity factors for electricity-districtheating IESs, which can estimate the effect of multiple coupling units on all power flows in the coupled DHS and EPS.They also discuss the benefits and limitations of PTDFs in IES.All methods proposing sensitivity factors for IESs, however, neglect the dynamic system behavior of IES as they derive the sensitivity factors from a steady-state power flow calculation.Hence, no method is available that represents the real system behavior.The methods do not allow a comprehensive analysis of the full potential of IESs (i.e., network storage and interactions) while having a high computational efficiency.
We extend the PTDF-based approach of [13] with a GS and enhance the approach by including the dynamic behavior of the DHS and GS.For this, we derive the sensitivity factors based on a joined quasi-steady-state power flow calculation presented in [14], which includes the dynamic behavior and network storage effects of the GS and DHS.With this, the proposed sensitivity factors include the full potential of IES.Thus, compared to PTDFs our sensitivity factors represent not only the power in the EPS but also temperatures, flow rates, and pressures in the DHS and GS.Our sensitivity factors allow a comprehensive investigation of possible effects on the overall system state when different units simultaneously adapt their power consumption or generation, including feedback effects between the different energy systems.Hence, no power additional power flow calculation needs to be run, reducing the computational burden.However, as the non-linear behavior of the system dynamics has a strong impact on the system state, we analyze if the sensitivity factors can represent the dynamic behavior of an IES.We present a detailed analysis in which we profoundly test the applicability of sensitivity factors for the use in IESs, their accuracy, and limitations.We test if sensitivity factors can be used to estimate the complex system behavior due to the system dynamics of DHSs and GSs and the interactions between the different energy systems.Finally, we compare the sensitivity factors derived for the quasi-steady-state power flow calculation with the sensitivity factors derived from a steady-state power flow calculation as presented in [13].
The main contributions of this article are as follows: r Deriving sensitivity factors for IES including the system dynamics of the DHS and GS r Analyzing the effect of the dynamic behavior of the DHS and GS on the sensitivity factors r Detailed investigation of the accuracy and limitations of sensitivity factors for IES r Comparing the accuracy of quasi-steady-state with steady- state sensitivity factors The article is structured as follows: Section II briefly describes the power flow calculation of the EPS, DHS, GS, and IES which is based on the Newton-Raphson method.Section III introduces the IES sensitivity factors while Section IV presents a detailed accuracy analysis of the new sensitivity factors compared to the results of a power flow calculation.

A. Electric Power System
As the EPS reaches a steady-state after a change within seconds the EPS is commonly represented as steady-state in the power flow calculation (e. g., [1], [15], [16], [17]).The state of an EPS is described by the complex nodal voltages, which are split into voltage angle δ N and voltage magnitude u N in the state vector x ps (e. g., [15], [16]): with the voltage magnitude u N being pu-values, based on the nominal voltage level of the investigated network U ref .
The nodal voltage angle and magnitude are determined based on the active and reactive power flow balances on all nodes except a slack node: Here, U N being a diagonal matrix of the complex nodal voltages u N .The node admittance matrix Y NN represents the network topology and the equipment characteristics.p p,N,set and p q,N,set are the pu-values of the known nodal active and reactive power injections or withdrawals (i.e., set points), which are based on a reference value P s,ref .
Both balances are joined in the vector of mismatches Δf eps .
A slack node must be defined with a known voltage magnitude and angle [9] as otherwise the Jacobian matrix is not invertible.

B. District Heating System
The DHS normally does not reach a steady-state within a simulation time increment of the power flow calculation due to the transfer delay of the temperature propagation.Thus, the dynamic behavior of the DHS must be included in the power flow calculation.The state of a DHS x hs is described by the mass flow rates on all edges in the system q m,E , the pressure difference of control elements Δπ CE (i.e., pumps at generators and differential pressure regulators at loads) and the temperatures at each node in the supply and return network ϑ N , which are affected by the dynamic thermal behavior of the DHS (e. g., [16], [18], [19]).The vector is of size (E + CE + N ) × 1, in which E is the number of edges, CE the number of pressure regulators, and N the number of nodes: In this, the vector of edge mass flow rates explicitly represents the supply and return network of the DHS.This allows a more realistic representation of DHSs as loads and suppliers are depicted as edges and the behavior of valves and pumps can be included [19].
A detailed description of the power flow calculation can be found in [19].In the following all variables describe the state in the current simulation time increment ν.If a variable depicts the state in a different time step the time step is given in the index.
The steady-state mass flow rates on all edges are determined by three equation systems (( 5), ( 6), ( 8)), resulting in equations equal to the number of edges E. First, the sum of all mass flow rates entering and leaving a node must be zero: in which E n is the set of pipelines connected to a node.The balance is set up for each node except a slack node.Second, the sum of pressure losses of each pipeline in each loop in the supply or return network must be zero: in which E m is the set of pipelines which are part of a loop.The balance does not include loops along loads or suppliers.The pressure loss along a pipeline is determined by: Here, L and D i are the pipeline's length and inner diameter, respectively, while Q m,l and ρ fl are the mass flow rate and the density of water.ξ is the pipeline's coefficient of friction, which is flow dependent and determined by the Poisson equation for laminar flows and the Colebrook-White equation for turbulent flows.Third, the heat flow rate at loads and supply units must be equal to a set heat flow rate Q th,d,set for each unit except a slack unit: in which c fl is the specific heat capacity of water while ϑ in and ϑ ex are the temperatures at the inlet and outlet of the load / supplier edge.As in the EPS, a slack generator is defined, balancing the DHS.
The pressure difference of pressure control elements (i.e., circulation pumps at suppliers and differential pressure regulators at loads) is determined by the simultaneous pressure control of [18].For this, a balance for each pressure regulator CE is set up as in (6) in which the sum of all pressure differences along edges which are part of the pressure control path must be zero.
The dynamic thermal behavior of the DHS affects the nodal temperatures as temperature changes propagate through the network.This propagation can be described by a simplified one-dimensional advection partial differential equation (PDE), Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.
which is of the form [20]: Solving the PDE to determine the outlet temperature of a pipeline leads to: L l (10) in which ϑ amb is the ambient temperature of the pipeline's surroundings and λ l is the heat loss per unit length of the pipeline.ϑ in,l,t−τ represents the temperature of the water element, that leaves the pipeline in the current simulation time step, when it entered the pipeline.Here τ is the transfer delay of the pipeline, i. e., the time it takes a water element to travel from the inlet to the outlet of a pipeline, which is defined as: In this, the transfer delay and the temperature loss depend on the dynamic behavior, i. e., the time a water element resides in the pipeline.The temperature propagation can be included into the power flow calculation by the temperature-gradient method described in [19], which shows a high accuracy up to a simulation time increment of 60 min.With this, an enthalpy flow rate balance is set up at each node, tracking the temperature distribution: in which Q h,n,in and Q h,n,ex are the enthalpy flow rates entering and leaving a node.E n,in and E n,ex define the set of edges entering or leaving a node.ϑ n represents the nodal temperature at the current time step.ϑ ex,i represents the temperature at the outlet of a pipeline and depends on the transfer delay and the temperature loss, i. e., the time a water element resides in the pipeline).
The above described balances are joined in the vector of mismatches Δf hs for the power flow calculation.

C. Gas System
The dynamic behavior of GSs is twofold and depending on the pressure level does not reach a steady-state within a simulation time increment of the power flow calculation.The gas compressibility has a strong effect in transmission systems while the hydrogen propagation affects both distribution and transmission systems.The state of a GS is described by the pressure at each node π N , the volume flow rate at each terminal of the network at standard conditions q v,n,Te , and the calorific value at each node h o,N [21].The vector is of size (2 N + T e) × 1, in which T e is the number of terminals: Here, the sign of the volume flow rate depicts if a flow is entering (positive sign) or leaving (negative sign) the pipeline.A detailed description of the power flow calculation can be found in [21].
The nodal pressures and terminal volume flow rates include the dynamic gas behavior arising from the gas compressibility, which is determined by six equations.First, the sum of all calorific value flow rates must be zero on all nodes, except at known pressure nodes, e. g., slack node: in which T e represents the set of terminals connected to a node while Q ho,n,dg,set is a known calorific value flow rate, representing a gas injection or withdrawal.The balance is set up as a calorific value flow rate balance to consider a varying calorific value due to hydrogen injection, ensuring that the heating demand at consumers is met.Second, a pressure balance is set up for all nodes at which the pressure level is known but the gas extraction or supply is unknown: in which π g,set is the set nodal pressure while π g,calc is the nodal pressure determined in the power flow calculation.Together both equations, ( 14) and ( 15), provide equations equal to the number of nodes N as each node has either a given gas demand/supply or known pressure.Third and fourth, a simplified continuity and momentum equation are used, adding equations equal to twice the number of pipelines L. The simplified equations assume a one-dimensional flow, compressible and homogeneous fluid, horizontal pipelines, and an isothermal flow.They are discretized by a fully implicit scheme in time and centered difference scheme in space.Choosing a space discretization equal to the length of the pipeline L leads to: Here, ν is the simulation time step while π l , Q v,n,l , and ξ are the mean pressure, the mean volume flow rate at standard conditions and the mean friction factor of the pipeline, which is flow dependent and determined as in the DHS.ρ n is the gas density at standard conditions while D i,l and A l are the pipeline's inner diameter and cross-sectional area.The isothermal speed of sound c is determined by the state equation: Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.
in which R and ϑ are the specific gas constant and gas temperature while the compressibility factor Z is determined by Papay's equation which is also used by [5] and [22].The equation is valid for pressures up to 150 bar and natural gas mixtures with up to 20 vol.-% of hydrogen [5].The gas properties of the gas mixture in each pipeline are determined by averaging the gas properties of hydrogen and natural gas according to their share in the gas mixture.
The fifth and sixth equation set are needed if the GS contains pressure regulators or compressors to determine the terminal volume flow rates.In this, we assume a steady-state volume flow rate through a compressor or pressure regulator: For compressors we also assume a dependency of the outlet pressure on the inlet pressure, i. e., the compression ratio, which is also done by [23]: In contrast, for pressure regulators we assume that the pressure is known at the outlet and that their gas demand is zero.With this we can apply (15).Equations ( 16)-( 20) add equations equal to the number of terminals T e in the GS.
As the nodal calorific values are affected by the dynamic behavior of the hydrogen injection, we use the calorific values to track the hydrogen propagation through the network, which is done in many studies (e. g. [5], [24]).The propagation of a change in hydrogen in a GS can be described by a simplified one-dimensional advection PDE [24]: which can be solved to determine the calorific value at the outlet of a pipeline as: As no source term is included in (21), the calorific value entering a pipeline reaches the outlet unchanged after a transfer delay τ along the pipeline: If the propagation is assumed to be steady-state (i.e., a change in calorific value occurs simultaneously on all nodes in the network) the transfer delay is zero.The dynamic behavior can be included into the power flow calculation by the calorific-valuegradient method described in [21].With this the calorific value flow rate balance is set up for all nodes except the slack node, tracking the hydrogen distribution: The difference between the nodal calorific value flow rate balance and the reduced demand/generation calorific value flow rate balance in (14) lies in the calculation of the calorific value flow rates.While ΔQ ho,n,DG,red assigns the incoming and leaving volume flow rates the same calorific value (i.e., the nodal calorific value), ΔQ ho,n,n assigns both flow rates different calorific values.The calorific value flow rate entering the node is determined by the calorific value at the end of the respective edge (first term in ( 24)), while the leaving calorific value flow rate is determined with the calorific value of the node (second term in (24)).The calorific value flow rate entering the node is determined by applying the calorific-value-gradient method described in [21] and thus considers the transfer delay of the hydrogen propagation.The calorific-value-gradient method shows a high accuracy up to a simulation time increment of 60 min [21].
The above described balances are joined in the vector of mismatches Δf gs for the power flow calculation.

D. Integrated Energy System
As the above-described power flow calculation methods are based on the Newton-Raphson method, the equation systems of the different energy systems can be easily joined, resulting in the state vector x ies and the vector of mismatches Δf ies : To improve the computational efficiency and to reduce convergence issues of the joined power flow calculation, the state variables in (25) are scaled to per-unit-values, reducing the order of magnitude of the values in the Jacobian matrix [25].In the EPS, the voltage magnitude u N is based on the nominal voltage level U ref .In the DHS, the nodal temperatures ϑ N are related to the minimum supply temperature of the generation unit at the slack node ϑ ref .Furthermore, the nodal pressures are based on the nominal pressure level of the network π ref,hs .In the GS, the nodal pressures π N and the nodal calorific values h o,N are related to the nominal pressure level of the network π ref,gs and the calorific value of natural gas H o,ref , respectively.
The Jacobian matrix J ies is set up based on the derivatives of the vector of mismatches Δf ies with respect to the state vector x ies : in which the submatrices on the main diagonal being the Jacobian matrices representing the single energy systems.The nondiagonal submatrices, on the other hand, represent the coupling and interdependencies between the different energy systems.For example, J h2p includes the effects of a change in the DHS on the EPS due to a change in heating set point of an electrode boiler (EB) while J p2g represents the effects of the EPS on the GS, resulting from a power set point change of an electrolyzer (ELZ).Our approach includes every coupling unit of the IES in Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.
the Jacobian matrix and not only units placed at the slack node of each network.Whether a non-diagonal submatrix contains non-zero elements depends on the coupling units in the IES and their operation mode.If the coupling units operate in heat-led mode, the heat generation is known while in power-led mode, the electricity demand/generation is known.The derivation of the submatrices are shown in more detail in [14].
Although the coupling units, such as heat pump, electrolyzer, electrode boiler and gas boiler (GB), have different working principles and connect different energy systems, their modeling can be generalized by a conversion factor f conv,cu : The generated power Q x,gen,cu can either be a thermal power Q th,cu , in the case of EBs and heat pumps, or a gas production Q ho,n,cu , in the case of an ELZ.The consumed power P x,con,cu , on the other hand, can either be an electric power consumption P p,cu , in the case of EBs and heat pumps, or a gas consumption Q ho,n,cu , in the case of a GB.Depending on the coupling unit and its operation mode (heat-led or power-led) either the consumption or generation must be known as otherwise the operation of the coupling unit cannot be obtained (e. g., [15], [16], [26]).
Although the conversion factor in ( 28) is shown as constant, the modeling approach also allows to include a power-dependent conversion factor.At each discrete time step ν of the power flow calculation the consumption or generation of a coupling unit is fixed and thus the conversion factor is constant in that point of time.However, between the discrete time steps, the conversion factor can vary depending on the power set point of the coupling unit.Such behavior can be included by using a different conversion factor in each discrete time step ν as either the power generation or consumption must be known for the power flow calculation and thus the conversion factor is also known.

III. SENSITIVITY FACTORS FOR IES
To estimate the effect of a unit's power change on the system state in an IES, we use PTDFs as they allow deriving easily the sensitivity factors based on the Jacobian matrix of the joined quasi-steady-state power flow calculation described in Section II-D.No further calculation or reformulation of the IES system state is needed as for other approaches, such as the Power Flow Decomposition [27] or Fractal Approach [28].
The PTDF approach, which is based on the Newton-Raphson method, linearizes the power flow balances at the system state of an energy system x ies,0 .Based on the linearization, the new system state after a change x ies,1 can be estimated: Here J −1 ies represents the sensitivity matrix, which is the inverted Jacobian matrix shown in (27).The Jacobian matrix can be derived from the Newton-Raphson method which is used to determine the initial system state x ies,0 .
To estimate the effect of a power change on the power flows in an IES, the power change is set in the corresponding element in Δf ies while all other elements are zero.Besides estimating the effect of a power change at any node in an IES, the sensitivity factors also allow the effect of a change in temperature in DHSs or calorific value in GSs to be analyzed.
Deriving the sensitivity factors from a joined power flow calculation allows the interaction between the different energy systems to be directly considered.In particular, if a change in one energy system affects the other energy systems, the feedback of the affected energy systems on the energy system in which the change occurs is included in the sensitivity factors.This would not be possible if the sensitivity factors would be derived from the independent power flow calculation of each energy system.A change in heat generation in the DHS, for example, not only causes a change of gas consumption and power generation of that unit but also a change of electric power consumed by circulation pumps.

IV. ACCURACY ANALYSIS
The following section presents two case studies investigating the accuracy of the sensitivity factors for different network topologies, placement of coupling units, as well as load and generation scenarios.In the first case study a city-district-sized distribution IES is modeled (Section V) while in the second case study a country-sized transmission IES is investigated (Section VI).In each case study different use cases are derived to analyze the accuracy of the sensitivity factor's estimation including a varying initial set point of the coupling unit, a variation in power change and different overall load situations.Each scenario is investigated for simulation time increments of 3 min, 15 min, and 60 min.The power flow calculation and sensitivity factor derivation is implemented in Matlab and the analyses are run on a standard office Laptop PC (Intel i5-7300HQ 2.5 GHz, 8 GB RAM).
To investigate the accuracy of the sensitivity factors, a power flow calculation is conducted for a given time horizon, e. g, 18 hours in the first case study and 42 hours in the second case study).Due to the varying energy demand and generation, a dynamic behavior in the DHS and GS arises over this time horizon.The last time step of the power flow calculation is used as the initial solution x pf,0 .Based on this initial solution, the sensitivity factors are derived.Then, only the generation of the coupling units is assumed to change in the next simulation time step while the consumption stays the same as in the initial solution.The effect of the coupling unit's power change on the power flow is estimated by the sensitivity factors and calculated by a power flow calculation x pf,1 .The deviation in % between the estimated change Δx sens and the actual change Δx pf is calculated by: V. CASE STUDY I: DISTRIBUTION SYSTEM The distribution system is based on the electricity-districtheating IES on Barry Island, United Kingdom which is widely Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.
1. Network topology of the distribution IES (DHS in red, EPS in green, and GS in yellow).The gray areas indicate the loads which are aggregated at the respective node in the EPS.Area 1 and 6 consist each of six MFHs, area 2 and 7 consist each of six businesses, area 3 consists of eight MFHs, area 4 consists of two MFHs and two businesses while area 5 consists of four MFHs and two businesses.used in the literature (e. g., [1], [16], [17], [26], [29], [30]).The network and load data ist taken from [14] (see Fig. 1).
The DHS consists of 20 loads representing multi-familyhouses (MFHs) and businesses, and has a temperature pairing of 110 • C / 70 • C (supply / return).Each load is located at a different node as shown in Fig. 1.The GS has a nominal pressure of 110 mbar and the same topology, load types and load profiles as the DHS.The EPS is represented as a 33 / 11 kV medium-voltage network.In contrast to the DHS and GS, the electrical loads are aggregated at the connection points to the low-voltage system (see gray boxes in Fig. 1).Furthermore, the distribution IES includes PV systems, a heat-led CHP unit, a heat-led GB, a power-led EB, a power-led ELZ and three circulation pumps, which are placed at the generation units in the DHS.The GB, EB, and ELZ have a conversion factor of 88 %, 99 %, and 60 %, respectively which are constant in the analysis.
The load and generation profiles are deduced by standard load profiles [31], [32] and the annual heating demand given in [30] (see Fig. 2).The rated heating and electricity demand can be found in [14].The EB and GB have a rated power of 300 kW each, reaching 18 % of the maximum heating demand in the DHS and GS.
The accuracy of the sensitivity factors is investigated by varying the power generation or consumption of the coupling units based on two use cases.First, the initial set point is fixed at 150 kW (50 % of their rated power) and the power change is varied between 0 kW and ± 150 kW.Second, the initial set point is varied between 0 kW and 300 kW while the power change is fixed to ± 10 % of the set point.In each use case, first the operation of only a single coupling unit is changed and then all coupling units are changed simultaneously.In the following analysis, the deviation between the estimated and actual change Fig. 2. Heating demand (left) [31], PV generation (middle) [32], and electricity demand (right) [32] of the distribution IES in pu.The red and green lines represent the profiles for the medium heating demand scenario while the gray lines represent the profiles for the high heating demand scenario.The electricity profile "G1" is the same in the medium and high heating demand scenario.
in system state are mainly presented for the high heating demand scenario.
In general, the following performance of the sensitivity factors was observed.Throughout all use cases, the estimation of the new system state using sensitivity factors is on average more than ten times faster than conducting a power flow calculation.The improved computational performance occurs as the Jacobian matrix only needs to be set up and inverted once to derive the sensitivity factors.In contrast, in a power flow calculation the Jacobian matrix is set up and inverted in each iteration.The accuracy of the sensitivity factors is generally better for the higher heating demand scenario and smaller simulation time increments.During a higher heating demand the power change of the coupling units have a smaller impact on the power flow.For smaller simulation time increments the dynamic behavior has a stronger impact on the system state as changes might not reach the end of lines.This mainly affects the nodal pressures, temperatures, and calorific values in the DHS and GS.Hence, the new system state is closer to the initial system state and the linearization error of the sensitivity factors is smaller.As the sensitivity factors estimate the new system state by linearizing around the initial system state, the estimated change has the same absolute value independently of a power increase or decrease.The accuracy analysis partly shows very large deviations of more than 1000 %, which appear when the absolute change of a state variable (e. g., nodal voltage, mass flow rate, calorific value) has a magnitude of 10 −4 or smaller.Such changes are mainly of numerical origins and would not affect system operation.Hence, only changes are considered in the accuracy analysis which have an absolute change above a threshold of 0.1 % of the reference value used for the power flow calculation.For example, a temperature change must exceed 0.8 K while a voltage change must exceed 11 V.The threshold is set according to the accuracy of measurement devices used in the different energy systems.

A. Power Change EB Only
If the EB's initial power set point is fixed at 150 kW and its power change is varied, the temperatures in the DHS change mostly on nodes connecting the main pipelines.In contrast, the temperature changes at loads are below the threshold of 0.8 K.The deviation between the estimated and the actual change in nodal temperature on the main nodes, however, is very large and increases linearly up to 50 % with an increasing power change ΔP (see Fig. 3, left).Also, the deviation increases the farther a node is located from the EB.The error of 100 % at a power change of 30 kW in the high heating demand scenario arises as the sensitivity is below the threshold while the actual change is closely above the threshold.Furthermore, in this case the accuracy of the sensitivity factors is slightly better for the medium heating demand profile which might be caused by the different system dynamics in the DHS.
The mass flow rates in the DHS have a similar behavior as the temperatures.While on pipelines connected to consumers a large error appears, the sensitivity factors provide a high accuracy for the main pipelines with errors mostly below 1 % for a simulation time increment of 15 min (see Fig. 3, right).
The changing heat generation of the EB leads to a changing operation of the CHP unit as a constant heat-to-power ratio is assumed.The sensitivity factors can determine the adapted operation of the CHP unit with a deviation below 2 % for the maximum ΔP and below 1 % in all other cases.The effect of the changed gas consumption of the CHP unit on the GS is estimated with a good accuracy which is mostly smaller than 5 % on pipeline 32, connecting the gas supply node with the CHP unit.For all other main pipelines, the error of the sensitivity factors is less than 22 %.The error decreases with larger ΔP as the change of volume flow rates increases up to maximum of 0.006 m 3 s and get farther away from the threshold of 0.001 m 3 s .The large errors at ΔP = ± 30 kW arise as the sensitivity factors estimate a change below the threshold while the actual change is slightly above the threshold.The deviation of the estimated and actual change in nodal pressures is quite large and lies mostly above 20 %.The absolute change, however, is small with approx.250 Pa being 2 % of the nominal pressure level.

B. Power Change GB Only
If the GB's initial power set point is fixed at 150 kW and its power change is varied, the estimated change of the sensitivity factors generally have a large deviation from the actual change.As the GB is located farthest away from the slack CHP unit, the GB has a large impact on the pipelines in its vicinity.Hence, the system state is much more sensitive to a change of the GB than to a change of the EB.In the high heating demand scenario, the deviation of the change in mass flow rates reaches up to 55 % compared to the actual change.This is a result from the small mass flow rate on pipeline 8 in the initial system state (see Fig. 4).Any change in mass flow rate on pipeline 8 leads to a new mass flow rate which is relatively far away from the initial solution.The farther away the estimate is from the initial system state, the larger the linearization error.Hence, the deviation in the medium heating demand lies mostly beneath 5 % for all main pipelines as the mass flow rate on pipeline 8 is larger.Nevertheless, the sensitivity factors can consider a mass flow reversal in contrast to [13] (see Fig. 4).
As in the EB use case, the GS is affected by a changed gas consumption of the CHP unit.However, due to the estimation errors of the mass flow rates in the high heating demand scenario no change in the operation of the CHP unit and thus no change in volume flow rates in the GS is estimated.In contrast, in the medium heating demand scenario the accuracy of the sensitivity factors is much higher as the estimation of the mass flow rate change is strongly improved.
The change of the GB's heat generation results in a changing gas consumption and thus changing volume flow rates in the GS on pipeline 30, 5, and 3 which connect the gas supply and the GB.For all other pipelines the change in volume flow rate lies beneath the threshold of 0.001 m 3 s .The sensitivity factors cannot predict the change in volume flow rate for a simulation time increment of 15 min and larger.The large deviation is a result of the estimation error of the mass flow rate which affects the estimation of the gas consumption of the GB.

C. Power Change ELZ Only
If the ELZ' initial power set point is fixed at 150 kW and its power change is varied, the DHS is not affected.In contrast to the EB and GB use case, the sensitivity factors can predict the change in calorific value at the ELZ (see Fig. 5, left).The estimation is independent of the calorific-value-gradient of the initial system state as a change in calorific value is directly introduced via the coupling unit and the nodal generation/demand calorific value flow rate balance.The deviation of the estimated and the actual absolute calorific value change is closest to zero at a power change of -90 kW as the actual change in calorific value shows a non-linear behavior (see Fig. 5, right).Due to the non-linear behavior which is a result of the gas compressibility, the difference between the estimation and actual change varies.The change in calorific value, however, is only predicted at the ELZ but not at the other nodes in the GS similar to the other two cases above.

D. Power Change All Coupling Units
Simultaneously changing the power generation or consumption of all coupling units results in a superposition of the results of the above described use cases.
If the set point of a coupling unit is varied and the power change is fixed, a similar deviation between the estimated and actual change state arises as in the above described analysis.
To show the effect of the dynamic behavior on the sensitivity factors' estimation we compare the proposed quasi-steadystate sensitivity factors with the steady-state sensitivity factors presented in [13].For this, the accuracy of the steady-state sensitivity factors are calculated by comparing their estimate with the quasi-steady-state power flow results.Comparing the accuracy of the quasi-steady-state and steady-state sensitivity factors shows that generally the quasi-steady-state sensitivity factors have a higher accuracy in estimating the system state after a change.The estimation of a change in mass flow rates in the DHS is up to 10 % more accurate with the quasi-steady-state sensitivity factors.Moreover, the accuracy in estimating changes in the volume flow rate in the GS is up to 80 % higher for small simulation time increments.With larger simulation time increments the accuracy of both sensitivity factor approaches become similar as the effect of the dynamic behavior decreases.Both approaches show the same accuracy in estimating the nodal voltages and nodal calorific values in all investigated use cases.This is because the EPS already reaches a steady-state for the investigated simulation time increments while the quasi-steadystate sensitivity factors are not able to represent the dynamic behavior of the hydrogen propagation.
The EPS has 11 loads and is represented as a 220 / 110 kV high-voltage network.The network is adapted by adding RES generation units in the 110 kV voltage level and additional nodes Fig. 6. 21-node transmission EPS used in the transmission IES, which is based on the IEEE-14 node test system presented in [33].Fig. 7. 18-node DHS used in the transmission IES, which is based on the 14-node DHS presented in [34].Fig. 8. 22-node Belgian gas transport network used in the transmission IES, which is based on the 20-node Belgian transport system presented in [35].
to include the coupling units of the IES.The DHS consists of 9 loads, each representing a city, and has a temperature pairing of 130 • C / 60 • C (supply/return).The GS consists of 9 loads representing different cities in Belgium and interconnection points to Luxemburg and France.The network has a nominal pressure of 55 bar and contains two motor-compressor stations, which have a compression ratio of 1.2 (CMP 1) and 1.5 (CMP 2) and an efficiency of 80 % [15].The network has two gas entry points of which node 1 is assumed to be the slack node.In contrast to [35], the GS only contains one storage at node 5.The transmission IES includes one heat-led and two power-led CHP units, a heat-led GB, a heat-led EB, two power-led ELZ, and five circulation pumps, which are placed at the generation units in the DHS.The location and the parameters of the CHP units, the GB, and EB are the same as in [15].The GB, EB, and ELZs have a conversion factor of 88 %, 99 %, Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.and 60 %, respectively which are constant in the analysis.The EB and both ELZs have a rated power of 25 MW.The EB and ELZs operate constantly at the rated power while for the CHP units the generation profile varies depending on the heating demand.
The load and generation profiles are deduced by the same standard load profiles as in the distribution IES (see Fig. 9) and the power given in [33], [34], [35].Additionally, PV systems are connected to nodes 10, 11, 12, and 14 which have a rated power between 0.5 MW and 2 MW while a wind farm is connected to node 3 replacing a fossil-fueled generation unit.
The accuracy of the sensitivity factors in the transmission IES is investigated as in the distribution IES.In the first use case, the initial set point of the coupling unit is fixed at 12.5 MW and the power change is varied between 0 MW and ± 12.5 MW.In the second use case, the power set point is varied between 0 MW and 25 MW while the power change is fixed by ± 10 % of the set point.Each use case investigates the impact of a change in power generation or consumption of the EB, both ELZs or all three coupling units.
The general behavior and accuracy of the sensitivity factors in the transmission IES is similar to the results of the distribution IES.

A. Power Change ELZs Only
If the initial power set point of both ELZs is fixed to 12.5 MW and their power change is varied, the sensitivity factors estimate a change in voltage magnitude with a high accuracy (see Fig. 10, left).Due to the linearization, the deviation increases with an increasing power change of both ELZs.Furthermore, the accuracy can be seen as independent of the simulation time increment as the behavior of the EPS reaches a steady-state for all used simulation time increments.
As in the distribution IES, the sensitivity factors only estimate a change in calorific value at the ELZ nodes.At ELZ 1, the accuracy is better for the high heating demand scenario as more gas is transported, and thus, the effect of a changing hydrogen infeed is smaller.At ELZ 2, the deviation for an increasing power change is approx.ten times higher than for a decreasing power change (see Fig. 11, left).These deviations arise from the non-linear behavior of the change in calorific value (see Fig. 11, right).As the non-linear behavior is less strong at ELZ 1 it is assumed that the non-linear behavior is a result of the gas compressibility, which has a greater effect on the pipelines around node 12.The observed deviation of 100 % at ELZ 2 at a power change of 0 MW is a result of the hydrogen propagation.Although no power change occurs, the hydrogen is still transported with the gas flow leading to a changing calorific value, which is not estimated by the sensitivity factors.
The gas compressibility also leads to a smoothing of the volume flow rate profile and to a delay of the change between the inlet and outlet of a pipeline, which is not represented by the sensitivity factors.Furthermore, the deviation strongly depends on the simulation time increment and the heating demand scenario as the effect of the gas compressibility is much more pronounced for large simulation time increments.

B. Power Change EB Only
If the initial power set point of the EB is fixed at 12.5 MW and its power change is varied, the sensitivity factors estimate a change in voltage magnitude with a high accuracy (see Fig. 10, right).Interestingly, the deviation does not have the expected V-shape but rather in a linear shape.Such behavior occurs as the sensitivity factors estimate a larger absolute change in nodal voltage magnitude than the actual change if the power is decreased while estimating a smaller absolute change than the actual change if the power is increased.In combination with (30) this leads to a larger deviation for a power decrease.The absolute difference between the estimated and actual change in nodal voltage magnitude, however, is the same for a decrease and increase.

C. Power Change of All Coupling Units
If the EB and ELZs are changed simultaneously, the deviation between the estimated and actual change of the system state has a similar accuracy as if the units are changed independently.The results can be seen as a superposition of the results of the above described use cases.
If the initial power set point of each coupling unit is varied and their power change is fixed, the overall results are similar to the above described use cases.
Comparing the results of the quasi-steady-state and steadystate sensitivity factors shows similar results as in the distribution IES.The accuracy of the quasi-steady-state sensitivity factors is generally better in estimating the impact on the mass flow and volume flow rates as well as on the nodal temperatures.In contrast, the accuracy for the nodal voltages and calorific values is similar.

VII. DISCUSSION
The sensitivity factors show a high computational efficiency in both case studies.Deriving a new system state is on average ten times faster compared to a power flow calculation.Besides their high computational efficiency, the sensitivity factors can accurately estimate the complex interactions in an IES and can estimate flow reversals.The sensitivity factors can estimate how power changes of multiple coupling units affect the operation of other generation units, in particular units connected to the slack node, and the effects on the power flows in the different energy systems.The accuracy, however, depends on the IES topology, current load and generation situation, and the unit's location and power change.Nevertheless, such a dependency is usual for PTDFs in EPSs due to the linearization of the system state [9].
The sensitivity factors can only estimate a change in system state if a power change of any unit in the IES occurs.Otherwise, the sensitivity factors assume the new system state to be equal to the initial system state (see (29)).Such behavior can be expected as the vector of changes Δf ies in (29) only contains zero elements if no power change occurs.The dynamic behavior of DHSs and GSs however, can lead to changing flow rates and thus a different system state even if no generation unit changes its power generation.Even if the power of a generation unit changes, it is difficult for the sensitivity factors to estimate the dynamic behavior of the DHS and GS.Due to the dynamic behavior, the IES responds differently to a unit's power increase and decrease (e. g., Fig. 11, right), resulting in a non-linear behavior of the IES.Even a small unit's power change can lead to a strong change in the IES' system state, which is far away from the initial system state.
Moreover, the gradient method used to include the dynamic behavior in the power flow calculation affects the sensitivity factors.If in the initial system state the temperature or calorific value gradients indicate an increase, the sensitivity factors will also estimate increase in temperature or calorific value.For example, if the previous calorific value gradient indicates a decreasing calorific value, the sensitivity factors extend this trend.This, in turn, will lead to a further decrease of the calorific value independent of the actual hydrogen injection at the estimated time step.
Interestingly, the sensitivity factors are not able to estimate the hydrogen propagation in the GS, although a temperature propagation in the DHS is estimated.Such behavior of the sensitivity factors is surprising as the hydrogen and temperature propagation are included using the same principles.An effect of the gas compressibility seems unreasonable as in the distribution IES a steady-state behavior was observed due to the small pipeline volume and pressures.A decoupling of the hydrogen propagation and volume flow rates might arise from the discretization of the gas flow PDEs, leading to the estimation error.This, however, could not be verified.Due to the above-mentioned restrictions the sensitivity factors are not able to fully represent the dynamic behavior of DHSs and GSs as well as the effect of the gas compressibility due to the non-linear behavior.Nevertheless, every approach deriving the sensitivity factors based on linearization will have problems predicting the non-linear, dynamic behavior.
In both case studies some use cases show a large deviation between the sensitivity factors' estimated change Δx sens and the actual change in system state Δx pf .The large deviations observed are a result of the calculation in (30) in combination with the actual change Δx pf being close to the threshold.In contrast, when comparing the estimated system state x sens and the actual system state x pf after a unit's power change, the deviation is less than 2 % in all use cases.The strong reduction in deviation is a result of the larger base value x pf , i. e., denominator in (30), and the small impact of a unit's power change on the overall system state of the IES.

VIII. CONCLUSION
This article extends Power Transfer Distribution Factors from electric power system analyses to integrated energy systems to estimate the effect of a power change on the power flows in the different energy systems.The sensitivity factors include the dynamic behavior of the district heating and gas system and consider the interactions between the different energy systems.For this, we derive the sensitivity factors from a joined quasisteady-state power flow calculation.We analyze the capability of the sensitivity factors to determine the effect of a power change on the system state and to represent the dynamic behavior of integrated energy systems using two case studies, comprising a larger distribution and transmission integrated energy system.The accuracy is determined by comparing the estimated change in the system state of the sensitivity factors to the actual change if the power generation and consumption of coupling units is changed.
The sensitivity factors show a high computational efficiency as they provide an estimate ten times faster on average than a power flow calculation while providing an accurate estimate Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.
of the interactions in an integrated energy system.The accuracy of the proposed sensitivity factors depends on the initial system state, the topology of the integrated energy system, the current load and generation situation, and the power change of the coupling unit.This dependency however, is normal for Power Transfer Distribution Factors.Although the sensitivity factors can provide good estimates, they are not able to fully represent the dynamic behavior of district heating and gas systems due to its non-linearity.Compared to steady-state sensitivity factors, however, our quasi-steady-state sensitivity factors have generally a higher accuracy as the dynamic effects are included in their derivation while having a similar computation time.
Considering the high computational efficiency, the highly complex interactions in IESs, and the dynamic behavior of the DHS and GS the proposed sensitivity factors are particularly well suited for network planning purposes in which many different scenarios must be tested.
The sensitivity factors' estimation might be improved by using the sensitivity factors in an iterative manner, similar to the power flow calculation but with less iterations or by including the history of gradients and nodal values in the Jacobian matrix.For both approaches, an analysis on the trade-off between the increase in computation time and increase in accuracy should be conducted.

Fig. 3 .
Fig. 3. Accuracy of the temperature at the EB (node 12) in the DHS estimated with the sensitivity factors compared to the actual temperature after the change for the medium and high heating demand scenario (left).Accuracy of the mass flow rate on pipeline 25 in the DHS for different simulation time increments and the high heating demand scenario (right).The accuracy is shown for different power changes of the EB and an initial set point of 150 kW.

Fig.
Fig. Comparison of the estimated and actual mass flow rates on pipeline 3, 8, and 10 in the DHS for the use case of a varying power change of the GB.The bar indicate the estimated mass flow rates while the lines represent the actual mass flow rates.

Fig. 5 .
Fig. 5. Deviation of the change (left) and absolute change (right) in calorific value at the ELZ (node 30) in the GS estimated with the sensitivity factors compared to the actual change in calorific value for the high heating demand scenario.The accuracy is shown for an initial set point of the ELZ of 150 kW, different simulation time increments, and different power changes of the ELZ (left).

Fig. 9 .
Fig.9.Heating demand (left), RES generation (middle) and electricity demand (right) of the transmission IES in pu.The colored lines represent the profiles for the medium heating demand case and the gray lines represent the profiles for the high heating demand scenario.

Fig. 10 .
Fig. 10.Accuracy of the estimated nodal voltage magnitude change at node 14 in the 110 kV EPS network for the use case of a varying power change of the ELZ (left) and EB (right).The accuracy is shown for different simulation time increments of the power flow calculation.

Fig. 11 .
Fig. 11.Accuracy of the estimated nodal calorific value change (left) and the absolute change (right) at the ELZ 2 (node 12) in the GS for the use case of a varying power change of the ELZs.The results are shown for the high heating demand scenario and for different simulation time increments of the power flow calculation while the absolute change is shown for a simulation time increment of Δt = 15 min.

behavior, integrated energy system, Power Transfer Distribution Factors, sensitivity factors.
q Flow rate in W, kg/s, kJ /m 3 , or m 3 /s.ϑ Temperature in • C. λ Heat loss in W /K m.