Modern Solutions to Modeling and Optimization of the Steady State of Multi-Energy Systems

. A complex solution of problem of creating mathematical models of multi-energy systems is possible using a unified approach. An approach that will ensure consistent construction of mathematical models and unification of computational algorithms. The paper presents the elements of the concept of energy circuits as a basis for unified modeling of systems of different physical nature. The existing and developed approaches to solving the problems of calculating the flow distribution in multi-energy systems are presented, based on the analysis of publications on this subject. The general structure of the mathematical model of the flow distribution in a multi-energy system and the set of optimization problems for steady-state operating conditions are described. A possible formulation of the optimization problem for the short-term operation of a multi-energy system is presented.


Introduction
Integrated management of an energy system means considering it as a multi-energy system taking into account the production, transmission, distribution and consumption of various energy resources in their interconnection. Management of the development (planning and design) and functioning (operation) of energy systems of this type implies the coordination of one-product (electrical, heat supply, gas transportation, etc.) subsystems included in them.
The main difficulties in the analysis and synthesis of multi-energy systems are associated with the heterogeneity of the observed physical phenomena and the complexity of their mathematical description and the variety of system elements. A complex solution of problem of creating mathematical models of multi-energy systems is possible using a unified approach. An approach that will ensure consistent construction of mathematical models and unification of computational algorithms. The basis for such a unified approach can be the concept of energy circuits [1,2], which in turn proceeds from the principle of the similarity of phenomena of different physical nature [3].
A mathematical model of the steady state of the system is a base of the calculation of the flow distribution. Such a model includes topological equations (corresponding to Kirchhoff's circuit laws) and component equations (determining the relationship between the state parameters along the arcs of the system graph) in its most general form for single-product system.
It is natural to apply a similar approach to multienergy (multi-product, i.e. with several energy resources) systems. However, when it is enough to take into account the laws of resource distribution and determine the values of the state parameters in case of singleproduct system, the calculation becomes more complicated in presence of several resources and their mutual transformation. The topologies of the system are different for different energy resources, and the interconnection (integration) of one-product systems is carried out through the nodes common to these systems -energy hubs [4]. Energy hubs are nodes that are converting energy resource from one type to another. Therefore, in addition to two or more types of energy resources, the model includes the equations for the interconnection of single-product systems (the equations of energy hubs).
In a multi-energy system (as well as in single-product systems) the calculation of flows distribution (determination of the values of the state parameters in the nodes and connections of the system that satisfy Kirchhoff circuit laws) has two purposes. First, it is performed to determine the existence and permissibility of operating conditions. Secondly, it is an essential component of solving the problems of states optimization in order to ensure the minimum cost of operation.
Traditionally, a means of such minimization is the redistribution of the load between sources, consumers and storage of energy resources. As a consequence, this redistribution applies to the elements of transmission and distribution networks. This means are supplemented also by the possibility of converting energy resources from one type to another in multi-energy system.
Covering the needs of the end consumer in different types of energy can be provided from different primary energy carriers through the chain (sequence) of their transfer and transformation from one type to another. Such a variety of the ways of energy supply, along with the possibilities of storing (accumulating) energy, provides necessary degrees of freedom for solving states optimization problems [5].
The paper presents the existing and developed approaches to solving the problems of calculating the flow distribution in multi-energy systems, basing on the analysis of publications on this subject. Besides approaches to formulations and methods of solving optimization problems of steady-state of multi-energy systems are presented.

Energy analogies principle and the concept of energy circuits
A characteristic feature of multi-energy systems is the heterogeneity of physical phenomena that should be described in the model. The concept of energy circuits [1,2] should be used in this case for a generalized representation of the components properties in technological chains of production and delivery of energy resources. This concept is based on the principle of similarity of phenomena of different physical nature [3].
Energy analogy means the identity of mathematical dependencies for systems of different physical nature. Considering this, a complex energy system can be imagined as consisting of a set of generalized components, which characteristics do not depend on the type of energy.
The existing theories of homogeneous in physical phenomena circuits (electrical, mechanical, hydraulic, pneumo-hydraulic, thermal, etc.) are the basis for the construction of mathematical models of the corresponding phenomena and technical systems. In its turn, the concept of energy circuits considers the general laws inherent to the energy circuit, independent of the properties of a particular type of energy resource. According to this concept, the conditions for the unification of the description of mathematical models of systems of different physical nature are as follows [1,2]:  Circuits must be formed from components defined using two types of generic parameters:  series parameters, which can be measured using a device series connected to a component and comply a law similar to Kirchhoff's first law,  parallel parameters, which can be measured using a device connected in parallel to a component and comply a law similar to Kirchhoff's second law.  The dimensions of both types of parameters must be physically compatible. It is preferable that the multiplication of the parameters determines the strength or flow rate of the energy resource.
In particular cases, generalized variables describe electromagnetic, mechanical, thermal and hydraulic phenomena, generating mathematically equivalent expressions. These expressions represent energy analogies (Table 1), due to the energy conservation law [6]. Possibilities of unification of the mathematical description and integration were studied in detail in [2]. This study considered different physical nature of elements, their properties, state and behavior. Also, this study considered the issues of analysis and synthesis of mechanism's drives s that structurally and functionally combine components of different physical nature, which require joint description (modeling). The concept of energy circuits was further developed in [6, 7, etc.].
For a circuit, represented by components of one physical nature, can take place internal transformations of energy from one form to another. These transformations are identical for electrical, mechanical and hydraulic circuits according to the principle of energy analogies. For each of these circuits, the energy conservation law can be represented as [2] ,  Adapted by [7] where and are parallel and series variables; is energy of parallel variable; is energy of series variable; ∫ is energy of dissipation.
There ⁄ is generalized capacitance characterizing the ability of a component to store energy in a form determined by the value of the parallel variable (kinetic energy of a moving body, energy of an electric field, potential energy of fluid pressure); ⁄ is generalized inductance, which characterizes a component's ability to store energy in a form determined by the magnitude of the series variable (potential energy of a solid, energy of a magnetic field, kinetic energy of a fluid flow). Thus, in accordance with the three terms of the energy conservation equation, resistive, capacitive and inductive parameters can be considered as inherent in energy circuits of any physical nature [2].
So, according to the principle of energy analogies, action variables determine the rate of energy dissipation in an element, and state variables are integral characteristics of action variables that change over time [7]. Then (Table 2):  rate of energy dissipation in an element is determined by the product of parallel and series action variables;  Kinetic energy is determined by integrating successive action variables over successive state variables;  Potential energy is determined by integrating parallel action variables over parallel state variables.

Model structure for calculating steady-state operating conditions
Traditionally, problems of calculating the steadystate conditions of the system include: 1) Calculating the flow distribution; 2) Determining the permissible steady state conditions.
The model includes technical restrictions of the controlled quantities. 3) Determining the optimal steady state conditions. The model contains an optimality criterion along with the flow distribution equations and constraints.

Flow distribution calculation problem
The mathematical formulation of the flow distribution calculation problem of a system is traditionally based on a topological model in the form of a structural graph. This makes it possible formalizing the process of obtaining a system of model equations in a generalized ordered form and using methods and algorithms of graph theory.
Hereinafter, we will consider the rate of movement of a resource in space and the rate of transformation of one type of resource into another as quantitative measures of the flow of energy resources. Respectively, the flow of energy resources is measured in units of mass, volume or energy per unit of time.
The mathematical model for the flow distribution calculation contains topological and component equations for each single-product system and the interconnection equations for common nodes of system. In addition, the model can contain the static characteristics of the injection along with potential changes ratios (this characteristic is used for transformer or throlling connections). The coefficients of the model equations are set as input data for solving the problem.
Topological equations (linear, corresponding to two Kirchhoff's circuit laws) include the constraint equations:  of series parameters in the node. It includes injections of the energy resource flow and medium flows in the incident edges;  of parallel parameters in the circuit. It includes active potential changes and passive potential changes in the edge resistances of circuit. Component equations (non-linear in the general case) connect the parallel and series parameters of each edge. It expresses the relationship between the flow in the edge and the passive change of the potential in its resistance.
The interrelation of models of single-product systems in it's simplest case is provided by the equations of system interconnection for nodes that are common to two or more systems (energy hubs). Note that, depending on the formulation of the problem, more detailed models of the hub can be used [4, 8-11, etc.].
In some cases, the process of calculating the flow distribution of a multi-energy system is advisable to be carried out alternately for single-product systems. Initial data (i.e. nodal injections of resource flow) of next calculation should be based on the results of previous one. Implementation of such an approach is presented in [12]. In this paper joint operation of electric and heat supply systems is being explored and the problem of calculating it is being solved.
When calculating the flow distribution, the number of unknown (dependent state parameters) is equal to the number of equations of the flow distribution model. The calculation consists in determining dependent parameters that satisfy these equations for given values of independent parameters It seems important to note the following consideration. In the general case, topological equations are a combination of balance equations corresponding to Kirchhoff's laws for nodes (the first law) and for circuits (the second law) of the system graph. In electric power systems, the application of the second law for calculating the flow distribution is associated with huge difficulties for most combinations of independent and dependent parameters of the state. This is due, firstly, to the complexity of the state parameters and, secondly, to the presence of power losses in the circuit elements. That is why for power flow calculations in electrical power systems Kirchhoff's first law is mainly used. If the parameters of the state are real and there are no losses during the transfer of energy resources (as in hydraulic calculations), then the two Kirchhoff laws can be used together or interchangeably. The obvious conclusion follows that nodal balance equations based on Kirchhoff's first law should be used for calculating the flow distribution as the most universal one.
Kirchhoff's circuit laws are linear. However, the mathematical model of nonlinear component equations (which is typical for most types of energy resources) will be a system of nonlinear algebraic equations.
There is still nothing better than classical Newton-Raphson method with some modifications for solving systems of nonlinear algebraic equations. As far as it can be judged from [13][14][15][16] and many other publications. Other methods proposed in the literature show greater efficiency in terms of convergence and speed. Usually, these methods have limitations on the scope themselves -either by the configuration of the graph schemes, or by the problems being solved. Therefore, they are inevitably combined with the Newton-Raphson method during algorithm formalization.
Such combinations seem to be justified because the Newton-Raphson method has a high sensitivity to the initial approximation. It means that, fast convergence of it takes place only with a "good" initial approximation. For example, the convergence of the simple iteration method is often far from the final solution. Therefore, a combination of several methods is an effective solution [17]. As an example: forward-backward algorithm [13,18] for calculating radial fragments of network, the method of holomorphic embedding of power flows [14] for finding a "good" initial approximation, optimization methods for obtaining fixed values of the state parameters (for example, the particle swarm method [15] is used for searching for "optimal" values of the active power of renewable sources). Thus, the Newton-Raphson method is currently the most effective method for calculating the flow distribution of various energy carriers systems, both single-and multi-energy, in combination with other methods that provide a "good" initial approximation. For example, in [16], this method was successfully applied to calculate the flow distribution of the combined gas transportation and electric power systems.

Determining permissible and optimal operating conditions
The calculation of the permissible conditions consists in determining the values of dependent and independent parameters. These parameters. These parameters must satisfy the following conditions: the equations of the flow distribution model are fulfilled, and the controlled values (dependent and independent parameters of the state, as well as functions from them) correspond to technical constraints in the form of equalities and inequalities. To be able to comply with the restrictions, the number of dependent parameters must exceed the number of flow distribution equations per number of restrictions.
It is necessary to change the independent parameters to ensure the existence of a solution if the permissible state does not exist at a certain step of iterative process.
State optimization problems are solved in conditions of operational and automatic states control of power systems and grids. The calculation of the optimal state consists in determining such an permissible state in which the objective function is equal to the minimum (or maximum) value. In this case, the flow distribution equations act as equality constraints on the state parameters. Costs for the time interval between two control actions, or costs per unit of time for a given consumer load at each moment of time are used as the target function. Variable cost components, which depends on the network conditions, are often taken as the target (minimized) function.
When optimizing the current state of the power sys-tem according to the power generation of energy resources, the volumes distribution of the produced resource per unit of time is determined. In this case, the producers' capacity injections that correspond to the minimum operating costs for given consumer loads are determined. A more particular problem is to optimize the electrical grid conditions. This means the calculation of the optimal node voltage, reactive power of sources and transformation ratios of variable transformers and autotransformers, taking into account technical constraints to minimize losses. In this case, the target function can be the losses of active power in the grid, which leads to the power flows distribution in the network with only active resistance taking into the account. Since the objective function is a quadratic form, and the constraints are a system of linear equations, in mathematical terms this is a quadratic programming problem. In the most general formulation, such a problem of optimizing the network conditions corresponds to the determination of the minimum active power of the balancing station and the damage to consumers from low-quality voltage [19]. The problems of optimizing the power system conditions and the grid conditions can be solved jointly (the so-called complex optimization of the operating conditions) [20]. However, separate optimization is used more often to simplify the task [21]: the power system conditions is optimized first, and then the power grid conditions is optimized.
The variables in all three described problems are continuous. Lagrange's method is usually used to account for equality constraints in the form of flow distribution equations. Non-linear programming methods are used to account for inequality constraints. Optimization of the combined cold-heat-power system is proposed in [22]. Reduction of the operation cost, impact on the environment, energy efficiency is used as the criteria. Three optimization strategies were used -electrical load follow, heat load follow and hybrid. Various optimization criteria (energy efficiency criterion, operating cost criterion, emission reduction criterion) were also applied to each of the strategies.
The model proposed in [11] calculates the amount of required energy carrier based on the current load and energy cost. In addition, it takes into account the presence of energy storage elements in the system.
In [7], many restrictions are described that are applied to the parameters of both the electrical system (capacitance of elements, nodal voltage, capacity of shunt capacitors, active and reactive power in nodes, active and reactive power flows), and water supply systems (pressure drop in lines and in units, step-down valves) and heat supply (water temperature, water consumption, heat transfer).

Energy hub modeling and optimization
The transmission and distribution infrastructures (grids) of a multi-energy system are different both in topology (which is obvious due to the different physical nature of the transmitted resource) and in the network configuration. These infrastructures are interconnected in common for them nodes (energy hubs), in which the transformation of energy types from one type to another takes place. Along with this, hubs are capable of producing, consuming and storing various types of energy carriers [5].
From a modeling point of view, the hub can be viewed as:  an interconnection point for models of systems of different physical nature; fixed nodal injection of resources in systems coupled through the hub (consumption of one type of energy carrier with the simultaneous generation of another of its types, taking into account the losses for transformation [7,9]);  an equation of the relationship between injections of different energy resources into the hub in the form of a static characteristic of nodal injections (terms are added to the static characteristic of the injection of a system node corresponding to a certain type of resource, which are functions of the node injection of another resource [5]);  complex system with its own flow distribution and capabilities to optimize functioning, and in which the conversion efficiency depends on the volume of the converted resource [10]. The authors of [10] propose an algorithm for calculating the flow distribution of an energy hub. In this algorithm, a preliminary calculation of all the parameters of the hub is done first, and then the data are transferred to the optimization unit.
In the case of multi-energy systems, when solving optimization problems, the authors most often pay attention to the connecting nodes of the systems, that is, to energy hubs. After finding all the parameters of the systems under study, the authors of [7] and [9] sought the optimal conditions, relying on injections into the energy hub. In the first case, it was a cogeneration turbine, and in the second, a gas turbine generator.
Optimization of energy hub control for short periods of time (20-30 hours) is considered in [11]. It is shows that the forecasting accuracy is mainly influenced by such parameters as changes in energy costs, fluctuations in demand, availability, and the effectiveness of the components. where -consumptions of the transported medium in gas, heat, cold and water supply; -nodal values of inflows and withdrawals of the transported medium.  Closing ratio -the law of variation of variables describing the functioning of the system at a moment in time : where -magnitude of change of variables. The following inequalities are taken as two-sided inequality constraints:  on power of generators,  on the bandwidth of the grid section,  to convert power from one type to another,  for the joint delivery of electric and heat power at the combined heat and power plant.
The solution of the formulated problem requires the development of new methods and development of existing methods. For this purpose, the existing approaches for calculating the flow distribution are used, as well as the procedures of mixed integer linear and nonlinear programming.

Conclusion
A characteristic feature of a multi-energy system is the heterogeneity of physical phenomena inherent in different types of energy resources. This creates difficulties in forming a unified mathematical model of such a system. A comprehensive solution to the problem of creating mathematical models of energy systems is possible within the framework of a unified approach. This approach ensures the consistent construction of mathematical models and computational algorithms.
The concept of energy circuits offers a generalized representation of the properties of the components of technological chains of production (generation) and delivery of energy resources to the consumer. Energy circuits are described using generalized series and parallel parameters. It is proposed to consider the rate of movement of a resource in space and the rate of transformation of one type of resource into another as quantitative measures of the flow of energy resources. In particular cases, such parameters describe electromagnetic, mechanical, thermal and hydraulic phenomena, generating mathematically equivalent expressions. This expressions also reflect energy analogies due to the fundamental law of conservation of energy.
The integration of single-product systems into multienergy systems is carried out through energy hubs, in which the type of energy resource is transformed. Such a transformation provides fundamentally new possibilities for optimizing modes specific to a multi-energy system. This important fact, apparently, explains the increased interest of researchers in the modeling and optimization of energy hubs and the immense growth of publications on this topic.
At the same time, the main purpose of the hub -the interface of two or more systems of production, transmission and distribution of different energy resourcesis often decreasing. The attitude towards research or optimization of hubs as isolated objects seems to be very controversial as seen in a number of publications. At the same time, the issues of modeling multi-energy systems in general have been worked out very poorly today and require deep and comprehensive research.
The equations of the flow-distribution model act as equality constraints on the parameters of the regime in optimization problems. It is possible to trace the continuity of the formulations and methods for solving problems of optimization of the regime of a multi-energy system with similar optimization problems for single-product systems. The main distinguishing feature of the optimization of the operating conditions of a multi-energy system is to determine the optimal combination of the types of energy resources used. So, the parameters of energy hubs are necessarily included in the number of optimized parameters.
The work was carried out within the framework of the State assignment under the Program of fundamental research of the Siberian Branch of the Russian Academy of Sciences, Project III.17.4.1 "Theoretical foundations of the creation and management of integrated intelligent energy systems" (reg. # AAAA-A17-117030310432-9).