Load-flow model of a multi-energy system

Current trends in the energy development includes the creation of multi-energy systems (MES) with several energy carriers that are planned, designed and operated with essential coordination of their subsystems. When managing the development and functioning of such systems, optimization problems should be solved to minimize the cost of production, transmission and distribution of energy resources. Traditionally, the means of such minimization are the redistribution of the loads of generation sources, consumers and energy storage devices. The topology of MES if different for different energy carriers, and topologies are combined through energy hubs – nodes in which energy is converted from one type to another. Such transformation can simply be considered as the consumption of one type of energy with the simultaneous generation of another type. This provides brand new opportunities for optimizing modes specific to a MES. Essential part of the search of optimal solution is the calculation of the load flow – that is, the determination of nodes’ and connections’ mode parameters of the system subject to balance equations (Kirchhoff’s circuit laws). The load flow equations act as constraints in optimization problems. This paper presents an attempt to form a single model for calculating the load flow of different energy


Introduction
In electric power, gas and heat supply systems, steadystate mode calculations (load flow calculations) in the network are carried out at all time and territorial levels of dispatch control, as well as in design and research organizations. Their purpose is to verify the admissibility of the planned regime in the short, long and long-term planning and solving a number of other issues. In addition, the methods and algorithms for calculating the steady-state mode of the electric network are the basis of methods for more complex calculations, in particular, methods for optimizing the network modes, comprehensive optimization of the power systems' mode by all variables, stability calculations, etc. [1].
Traditionally, energy (electric, gas, heat supply) systems are designed and operated separately. Meanwhile, the interdependence between these systems is becoming increasingly apparent, for example, power generation is increasingly dependent on the gas supply system. In addition, the composition of the problems of managing the operating modes is in many respects similar for the production / transmission / distribution / consumption systems of one or another energy resource. In particular, when managing the modes of electric power, gas and heat supply systems, such problems are solved as state estimation (forming a model of the current mode), steady state calculation, operational adjustment of optimal operating modes within an acceptable area, analysis of possible emergency situations, a comprehensive assessment of the reliability of networks, operational control of parameters and admissibility of the current regime, assessment of reliability indicators of the current mode, assessment of mode's technical and economic indicators, the adjustment of modes and networks by the conditions of reliability, the restoration work plan formation after accidents etc.
Note that all of these problems require load flow equations to be calculated. The mathematical formulation of the problem is concretized at the stage of clarifying the temporal level of mode control taking into account external conditions, the actual state of the system, the availability of the necessary initial data, and other factors.
One of the possible ways to solve the problem of rational distribution of energy resources is the creation of multi-energy systems. The purpose of such systems is the combined production, transmission, conversion and storage of several types of energy resources. Accordingly, a whole range of control and coordination problems is associated with this, both at design and operation stages.
Currently, the joint control of the modes of generation of electric and thermal energy is being held at the level of cogeneration sources -cogeneration plants. At the same time, the issues of controlling the regimes of trigeneration and, especially, quadrogeneration plants need serious theoretical and methodological study. To an even greater extent, this applies to the management of the functioning modes of MES. Meanwhile, a joint analysis of, for example, gas and electric infrastructure would make it possible to evaluate their productivity, the economically optimal load on the gas and electricity supply network, satisfying the physical restrictions on the regime parameters while satisfying demand [2].
The load flow problems of various energy resources are also similar from the calculating point of view. Quantitative characteristics of the energy carrier (substance) flow or energy are determined in all problems, which is the speed of their movement in space or the speed of their conversion from one type to another. The substance flow (more exactly, the flow rate) is measured in units of mass or volume per unit of time, and the energy flow (electric power) in units of energy per unit of time, respectively. Systems of nonlinear equations are compiled in all problems in one of two (interchangeable) forms -in the form of a balance of flows in nodes or in the form of a balance of potential changes in the contours of a scheme. The mathematical methods being used in this case are also almost identical.
The difference caused by the different physical nature of the considered energy resources can be traced only in the dependences between the flow along the edge of the graph of the design scheme and the change in potential in the resistance of this edge. These dependencies generally are nonlinear and that is reason why the system being solved is nonlinear.

Optimization problems in a multienergy system
The purpose of applying methods for calculating optimal modes is primarily to find acceptable, i.e. satisfying the conditions of reliability of power supply and proper quality of energy and modes. Of no less importance is the choice of the most profitable modes, since this allows ensuring savings by reducing fuel consumption and network losses for practically no additional costs [1].
Satisfying the needs of the end consumer in different types of energy can be ensured from different primary energy carriers through a chain (sequence) of their transfer and conversion from one type to another. Such variability in the ways of providing energy supply to consumers along with the possibilities of storing energy grants the opportunities for solving problems of energy supply optimization.
The problem of estimating the optimal mode is to find the optimal values of all parameters characterizing its affordability and efficiency. Comprehensive optimization of the regime is required for all variables, which involves the selection of optimal values of the controlled parameters, which, when optimized, are considered as independent variables [1].
The mode optimality criteria is the minimum of total energy costs (or total operating costs) while meeting the requirements of energy and energy carriers demand and observing all technological restrictions for a given period of time [2,3].
Variable constraints include: • Equality constraints corresponding to the initial data of the load flow problem; • Equality constraints that coincide with the equations of the mathematical model of the load flow problem; • Bilateral inequalities as the results of solving the load flow problem. From a math point of view, the problem of multienergy optimal load flow is considered as a nonconvex nonlinear bounded optimization problem.

Math model for load flow calculating
The math formulation of the load flow problem is based on a topological model in the graph form, which makes it possible to formalize the process of obtaining the set of equations of the model in a generalized ordered form. This representation allows the use of methods and algorithms of graph theory [4].
In load flow calculations, energy sources can be associated with both nodes and edges of a system graph. In the first case, the sources are represented by nodal injections, and balance equations are compiled in accordance with Kirchhoff's first law. In the second case, the sources appear to be active potential difference and the balance equations of the Kirchhoff 's second law are used. In particular, in hydraulic load flow problems, nodal injections and "active" pressure differences at edges are traditionally taken as input, and resource flows and "passive" pressure differences at edges are considered as variable (calculated) values [5].
The structure of the studied system is described by a directed graph containing N nodes and L edges. The incidence matrix (connection matrix) of this digraph can be expressed as (1) Elements of this matrix are equal to one if the edge l leaves the node n, -1 -if it enters into it, 0 -if the edge is not incident to this node. The incidence matrix fully describes the configuration of the graph.
Each l edge of the digraph is associated with its resistance (system parameter) . Here and below n and k are nodes incident to l edge and double lower index nk means the orientation of the edge from the node n to the node k. Upper and lower indices will be skipped further in case of insignificance for the current study.
The vector of edges' resistance is Together with the matrix , this vector completely determines the topology of the system (that is, along with the configuration, also the "weight" of the edges).
It is necessary in the traditional method of load flow calculation to set a number of mode's parameters. The problem comes down then to determining the remaining mode's parameters in nodes and edges. The load flow equations are considered as constraints in optimization problems. Since optimization requires one or more degrees of freedom, this means that one or more parameters that were fixed in the load flow problem are "released".
The system mode parameters matched to each node n of the digraph: • Injection of the transported resource flow , which matches with the nodal flow rate of the substance (volume or mass of the transported liquid or gas per unit time) or injection of electric current (or electric power). Injection is considered to be positive when importing a resource into the system and negative when exporting it from the system. Thus, a positive injection matches to the generation of a resource, and a negative injection corresponds to its consumption; • Nodal potential difference , corresponding to nodal pressure (or its square) in hydraulic networks and nodal voltage -in electric ones. • The system mode parameters matched to each edge l of the digraph: • The flow of transported resource , which matches with the flow rate of the substance in hydraulic networks, and with current or power in electric networks. According with the Kirchhoff's law, the sum of the flows along the edges incident to the node n is equal to the injection of this node: ; • "Passive" potential difference , which matches with the pressure difference or a voltage drop in the edge's resistance: ; • "Active" potential difference due to the potential source on the edge (which matches to the pressure difference on the edge during the operation of pumping units or electromotive force).
All these parameters (resistance of edges, nodal injections and potentials, flows and potential differences on edges), in the general case, can be considered as complex quantities.
The vectors of nodal injections and potentials for nodes can be expressed as: and vectors of flows, passive and active changes in potential -for edges: Taking into account the introduced notation, the mathematical formulation of the load flow problem can be described as (5) (6) (7) Equation (5) expresses the flows' balance (the Kirchhoff's law) at all nodes. Equation (6) determines the balance between the potential difference at the final nodes ( ) and the algebraic sum of the passive and active potential differences ( and ) of arcs. Equation (7) expresses the dependence between the flow and the passive potential difference in the edge's resistance l (such dependencies should be determined for all ). Equation (7) is nonlinear in the general case -with the exception of electrical calculations (where the flow matches with electric current), as well as hydraulic calculations during laminar motion of the substance. In these two cases, the dependence between the flow and the passive potential difference in the edge's resistance is linear and can be expressed as: Quadratic dependence between the friction pressure loss and edge flow (turbulent motion of the substance) is the most common dependence in hydraulic calculations of piping systems through which an incompressible medium is transported and can be described as: The passive potential difference is an increasing continuous flow function, taking values from -∞ to + ∞ and equal to zero at zero flow.
In electric power calculations (when the resource flow matches with complex electric power) equation (7) takes the form: The dot above the letter symbol indicates the complexity of the indicated quantity, and the sign indicates the conjugation of the complex quantity.

Energy hub model as a static characteristic of nodal injection
Network models of different energy carriers are interfaced through energy hubs, which are energy converters of networks. The energy hub is the link between production, transmission and distribution systems. It exchanges energy with sources, consumers, and networks through hybrid input and output ports. A typical energy hub consumes various types of energy at its input ports (connected, for example, to electric, gas supply and heating networks), and produces some types of energy (such as thermal and electric) at the output ports. In electric load flow calculations, nodal injections are often represented by functions of variable mode parameters. For example, the static characteristic of the of node n injection as a function of its voltage can be represented as a quadratic dependence with constant set coefficients (some of them may be equal zero): The terms on the right-hand side represent the power of a fixed-conductance shunt, a fixed current source, and a fixed-power source (shunt) respectively.
A model of an energy hub can be formed in the same way. The transformation of the type of energy carrier (converting gas to steam or electricity, steam to electricity, diesel fuel to electricity, electricity to heat, etc.) can be simplified as the consumption of one type of energy carrier with the simultaneous generation of another type (taking into account conversion losses). The terms that are the functions of the nodal injection of another resource are added then to the static characteristic of the injection of the system node corresponding to a certain type of resource.
For example, if the node n of the system is a hub containing a gas turbine unit (GTU), then for a positive nodal injection in the electrical part of the hub (electric power generation) addition ( ) can be represented as a function of negative injection (resource consumption) in its gas part: is an injection at the node n of the gas network, is a specific gas consumption for the production of a unit of electricity, is a fraction of injection used to generate electricity. At the same time, the term is added to the static characteristic of the injection of node n of the gas network and can be expressed as: There are may be more terms in the static characteristics, since the total injection of the electric part of the hub can consist of both the above generation and consumption (for example, for the thermal energy production).
The fraction of the resource injection that is converted to another type of resource indicated in the equations presented here will play the role of an adjustable parameter (independent variable) in optimization problems of the multi-energy system. So, for a gas network, gas turbine is considered as a negative injection (gas consumer), the value of which depends on the mode of the electric network. On the contrary, gas turbine units seem to be a positive injection (source of electric power) in the electric network.
The simplest model of energy or energy resource storage is a positive or negative injection, depending on its operating condition (output or consumption of a resource).

Regarding load flow calculation methods
Sources (resource flow and active potential change) can be located both in nodes and on the edges of the graph in equivalent scheme of a multi-energy system, as can be seen from Section 3. This arrangement determines the form of the balance equations of the mathematical model. At the same time, the placement of potential sources in nodes can lead to an unjustified increase in the number of nodes of the design scheme.
On the other hand, since in this paper we use the concept of energy hubs (as common nodes for networks with different energy carriers), we have to put up with this. For example, if we consider a gas pumping unit (GPU) not only as a source of active pressure, but also as a consumer (for example, electricity), then it will have to be placed in the graph node of the gas supply scheme. However, if the GPU is powered by a local source, it is advisable to place it on the corresponding edge of the graph.
Assuming that in equivalent schemes both methods of source arrangement would have to be used, the method of load flow calculation should be combined. This method combines different locations of sources in a single mathematical model [5]. The method is based on the representation of the vectors of the regime parameters and the incidence matrix in block form: where the subscripts D and V represents, respectively, the given and sought quantities. The equations (5) These four matrix equations, along with a set of coupling equations (7), represent a load flow model for an equivalent scheme with a part of the sources are located on edges and the other part -on nodes.

Conclusions
The function of the energy system (regardless of which energy carrier is implied) is the power supply of consumers. A multi-energy system is a system with different and, at the same time, interchangeable energy carriers. So, the task can be set to find a rational or optimal (by some criteria) combination of energy carriers used.
The introduction of multi-energy systems will improve the whole range of performance indicators of the energy system, the main of which are optimization of consumer load schedules, optimization of energy production modes and network infrastructure, improving energy quality and reliability of energy supply.
The infrastructures (networks) of transmission and distribution of a multi-energy system are different both in topology (which is obvious due to the different physical nature of the transmitted resource), and network configuration. These infrastructures are interfaced between themselves in nodes common to them (energy hubs), where the transformation of types of energy from one type to another takes place. Along with this, hubs are capable of producing, consuming and storing various types of energy carriers.
Depending on the specific optimization problems, the load flow calculations of different energy carriers can be integrated at the levels of: • mathematical model (system of nonlinear equations), • steps (iterations) of each calculation, • calculation results (in the general case, iteratively).
The load flow calculation of energy carriers (the calculation of the steady state based on Kirchhoff's laws) is the basic one for any statement of the problem. Methods, algorithms and software for such calculations are developed for specific energy carriers, but still not developed for combination of them.
This paper focuses on the first of these coupling levels. A complex load flow model in gas, heat and electricity systems is proposed. The model is a system of nonlinear algebraic equations, including flow balance equations at all nodes of the structural graph of the energy system (corresponding to the Kirchhoff's law), balance equations of the potential difference at the end nodes and potential changes on the edges, as well as the equations of the connections between the flow of the medium (or energy) by the graph edge and the potential difference of the vertices incident to it.
It is proposed to simulate energy hubs in the load flow calculations of MES in a simplified way of sets of interdependent static characteristics of parts of the hub related to different energy resources.