Aggregation of distribution heating networks in the problems of hierarchical optimization of large heat supply systems

The article shows that the problem of dimension is a significant obstacle to the automation of solving the problem of finding optimal hydraulic modes of heat supply systems. To overcome this problem, the authors have previously proposed a hierarchical approach to optimizing the hydraulic modes of heat supply systems. The article discusses the problem of aggregating distribution heating networks within the framework of this approach. To solve it, the previously proposed by the authors of loop reducing dynamic programming method is adapted. Its operability and computational efficiency are checked on computational experiments.


Introduction
Energy efficiency problems are topical both in Russia and abroad, while heat supply systems (HSS) have significant reserves of energy saving [1], which can be released by organizing optimal modes of their operation. In practice, the task of planning HSS modes is solved by multivariate calculations of the mode [2]. The choice of methods for organizing the modes is assigned to the specialist performing the calculations, which does not guarantee the optimality of the solutions obtained. The automation of solving these problems is complicated by a number of factors: the high dimensionality of the HSS [3], [4], the nonlinearity of the involved flow distribution models, the presence of several objective functions, etc. For these reasons, there are no methods and software systems suitable for wide practical application. This determines the relevance of the development of methods for calculating the optimal HSS modes.
One of the main problems arising in the optimization of HSS modes in large cities is associated with their large dimension, calculated in many hundreds of thousands of nodes of the design scheme. In a number of works, for example, [5], this problem is not solved, limiting itself to considering only low-dimensional HSS. In some works ( [6,7] and others), an approximation of the dependence of the value of the objective function on the parameters of the mode selected as a basis is used, which significantly complicates the correct accounting of discrete variables associated with the composition of the equipment. In other cases, to overcome the dimensionality problem, aggregation of HSS schemes is used, which prevents taking into account the entire set of constraints and does not guarantee the required accuracy of solutions. A common approach is using ready-made solvers. The main disadvantage of this approach is the impossibility of adapting methods to specific problems, and, as a consequence, too high computational costs. So, for example, in [8,9] the physical model of the network is built in the Simulink / Matlab environment, for calculating the state of which the CPLEX solver is used, and for optimization -the ReMIND software. The use of genetic and evolutionary algorithms requires even greater computational costs. For example, in [7], a nested iterative cycle is used to minimize the cost of pumping the coolant. In the inner loop, the allowed mode is calculated using the SIMPLE algorithm. On the outer loop, the parameters determining the mode are changed using a genetic algorithm.
Mainly, problems related to operational management are considered, for example, [10 -12]. Thus, at present, there is no satisfactory solution to the problems of planning HSS modes of real dimensions, which arises at the stage of preparation for the heating season.
The subject of this article is the tasks and methods of optimizing the operating modes of the HSS. The object of application is distribution heating networks (branched HSS fragments to end users) as part of hydraulically connected HSS in the process of optimizing HSS modes. It is assumed that the temperature conditions at heat sources (HS) are set, the heat losses in the networks are eliminated, and their residual value can be neglected. In this case, the requirements for providing consumers with thermal energy is reduced to the need to maintain their required coolant flow rates, and the task is reduced to the optimization of the hydraulic mode (HM).

Statement of the problem
Meaningfully, to optimize the HM of the HSS, it is necessary to find control actions that implement the mode that meets the admissibility requirements and achieves the specified optimization goals.
Earlier [13], the authors proposed a hierarchical approach to optimizing HSS modes to solve the dimension problem. This approach includes the following stages: 1) decomposition of HSS into backbone (BHN) and distribution (DHN) heating networks; 2) search for the limits of permissible change in the mode parameters at the decomposition point; 3) optimization of the BHN mode, taking into account the constraints obtained in the previous step; 4) optimization of DHN modes taking into account the parameters of the BHN mode at the decomposition point. BHN contains all heat sources, pumping stations and a single-line multi-loop part of the network. DHN includes only branched passive networks to end users. The decomposition point is the point of joining the DHN to the BHN in a single-line representation and two nodes in the bilinear representation, one of which (point a) is the junction of the supply pipelines BHN and DHN, the second (point b) -the return ones ( Fig. 1). For BHN, the decomposition point is a generalized consumer with bilateral restrictions on pressures in the supply and return pipelines, as well as on the difference in these pressures. For DHN, the decomposition point is a generalized HS modeled by two dangling nodes. Thus, when aggregating the DHN (replacing the DHN with a generalized consumer), it is necessary, based on the information about the DHN, to determine the minimum and maximum pressures at points a and b, as well as the minimum and maximum difference between these pressures at which at least one DHN hydraulic mode that satisfies the admissibility requirements.

Model of controlled flow distribution.
For DHN, due to the branched configuration, there is a fixed, easily computable flow distribution. Accordingly, the model of controlled flow distribution [14] will take the form: where: A -is the matrix of incidents of the DHN design diagram; P -m-dimensional pressure vectors; x, y -n-dimensional vectors of flow rates and pressure drops on the branches; h(x,z) -n-dimensional vector function with elements hi(xi,zi), i I  , approximating the dependence of the pressure drop on the flow rate on the branches; z is the vector of relative hydraulic resistances (for example, increased by throttling); X=(P, x, y, z). As an approximation of the hydraulic characteristics of the i-th branch, we can take [15]: Here: si -hydraulic resistance. In (2), to simulate a control action mode that is absent or forbidden for change at the planning stage, the corresponding variable can be equated to a constant.
If any mode parameter depends on the external environment, the corresponding variable is fixed. Let us call these parameters boundary conditions. For HSS, such parameters are the pressure in the make-up nodes (the nodes in which the coolant is supplied to the HSS to compensate for coolant losses and withdrawals) and the flow rates in all other nodes.
Let us write down the requirements for the admissibility and feasibility of the mode as [14] . Then the system of restrictions takes the form: When solving the problem of DHN aggregation within the framework of the hierarchical approach, the following are considered known for each DHN: topology of the calculation scheme; border conditions; coefficients of hydraulic characteristics (si, i I  ); permissible limits of variation of continuous variables. You need to find the following quantities: min a P under constraints (1) and (3);  (1) and (3).
Here a is the junction of the supply pipelines BHN and DHN, and b is the return ones.

Loop reducing dynamic programming method
Earlier [15], the authors proposed a loop reducing dynamic programming method (LRDP) to solve the problem of optimizing the DHN HM. The advantages of this method include: guaranteed finding the optimal HM; better performance in comparison with alternative methods (for example, based on the method of interior points [16]); growth of computational costs linear in the dimension of the problem; the possibility of solving the problem of multiobjective optimization in one application of the method. The main idea of LRDP is to consider the whole DHN. To facilitate tracking the second Kirchhoff's law compliance along the contours during the forward stroke, the contours are reduced to the state of two branches connected in parallel during the reduction of the entire design diagram to one branch with cutting off both unacceptable and non-optimal pressure distribution options. Instead of the recurrence relations typical of the dynamic programming method used for branched pipeline systems [17,18], to preserve the possibility to apply the principles of dynamic programming on "reduced" loops, special techniques are used for equivalent sequential and parallel connection of branches, aggregating pressure drop trajectories. On the return stroke, the only remaining pressure distribution option is restored to the entire network.
The intervals of the allowable change in the pressure are divided into equal nonintersecting subintervals (pockets). Let's renumber the pockets so that their numbers represent discrete pressure readings at the nodes with some factor. For the sake of definiteness, assume that all branches are directed downstream. Let's designate the initial node of the i-th branch as g . This representation of the DHN branches allows us to consider a fragment consisting of several sequentially connected branches as one branch. You can also consider as a branch a fragment of an DHN, consisting of several parallel-connected branches.
The direct stroke of the LRDP consists in contracting the trajectories of the pressure drop with rejection of unacceptable and non-optimal variants on the design diagram of the DHN due to the equivalenting techniques, as a result of which there will be only one branch with the optimal pressure distribution.
To equate parallel branches, do the following. If on two branches connected in parallel (i1, i2), there are a pair of segments (one per branch), which coincide with the initial and final pockets, respectively ( 1 ), then such a pair becomes a segment on the equivalent branch, starting and ending in the corresponding pockets. The value of the objective function increment for it is equal to the sum of the objective function increments of both equivalent segments. All segments that are not included in any such pair are discarded.
When equating a fragment of an DHN, consisting of two branches (i1, i2) connected in series (for definiteness, fi2 = li1), with one branch, the following is done. There are all such pairs of segments on these branches ( 1 2 1 2 , k k i i g g ) that at the node fi2 these segments have a common pocket ( 1 ). If two pairs of segments are found connecting the same pockets, the pair with the worst value of the objective function increment is discarded. If after rejection there are several pairs left, any of them is selected. Each of the found pairs of segments on the equivalent DHN fragment turns into a segment on the equivalent branch with the corresponding initial and final pockets. The value of the objective function increment for it is equal to the sum of the objective function increments of the segments of the equivalent pair. The reverse stroke consists in restoring the optimal trajectory of the pressure drop (vector P) along all the elements of the original design diagram of the DHN. According to the known vector P, as mentioned above, it is possible to restore the remaining parameters of the optimal DHN HM. To restore the optimal trajectory of pressure drop, it is necessary to remember for each equivalent branch which fragment it is equivalent to, and for each equivalent segment -which segments it equivalent. Knowing this information, the segment remaining after the DHN design diagram reduction into one branch can easily be expanded to the entire DHN design diagram.

Search for the limits of changing the mode parameters at the decomposition points
At given pressures at points a and b, it is possible to search for an acceptable mode using LRDP, taking 0, , In this case, when equivalenting branches connected in series, instead of choosing the best pair of segments for two pockets, one can choose any pair for these pockets that satisfies the condition . This reduces computational costs.
When searching for the optimal or permissible HM of the DHN, at points a and b, the pressures caused by the HM of the BHN are set. When searching for the limits of the permissible change in the parameters of the mode at the decomposition point (points a and b), at these points the pressure change intervals are set due to technological limitations.
To solve the problem of DHN aggregation, it is proposed to apply the following computational scheme of the LRDP.
1. For all branches of the DHN, all allowed (within the constraints on the admissibility of the mode) segments are determined.
2. If on all branches the sets of allowed segments are not empty, go to step 4. Otherwise, there is no solution, exit.
3. Find all fragments of the DHN, consisting of sequentially connected branches. Each of them needs to be equivalent to one branch using the described technique for equivalent sequential connection of branches.
4. Find all fragments of the DHN, consisting of parallel-connected branches. Each of them needs to be equivalent to one branch using the described technique for equivalent parallel branches.
5. If the DHN is reduced to one branch, go to point 7, otherwise to point 3.
6. There will be the set of segments on the remaining branch. In this case, for each pair of pockets at points a and b, the following condition will be met: if there is at least one admissible HM for a given pair of pressures, in the found set there is a segment corresponding to the admissible HM and connecting these pockets. It remains to search for the minima and maxima of the explored parameters for this set.
The main advantages of this method are: no iterative processes; no need to solve large systems of equations; the ability to solve problems (4) -(9) simultaneously in one pass.

Computational experiment
The purpose of the computational experiments was: 1) checking the possibility of the proposed method to find solutions to problems (4) -(9) simultaneously; 2) comparison of the performance of the proposed method in solving problems (4) - (9) with the performance of continuous optimization methods [13].
As reference methods for solving problems (4) -(9), we used continuous optimization methods [13], which are guaranteed to find solutions to these problems. The essence of these methods lies in the iterative reduction of the uncertainty interval, which obviously contains the optimal values of the objective function (the investigated quantity). At each iteration, the search for an admissible HM was performed by the method of interior points [19,20], which is an iterative process, at each iteration of which a system of algebraic equations of large dimension is solved.
To test the operability and efficiency of the proposed method, a series of computational experiments were carried out. Below is a typical example of a computational experiment.
In the computational experiment, problems (4) -(9) were solved on a conditional DHN, the diagrams of which are shown in Fig.2   For the case with the same restrictions, the following restrictions were obtained: 1 74 P  , 1 120 P  , 16  The verification of the correctness of the solution of problems (4) -(9) was carried out by reference methods. The check showed the correctness of the results of the proposed method. The time spent on the search for the limits of changing the mode parameters at the decomposition point by the proposed method (less than 5 minutes) was less than the time required to solve this problem by continuous optimization methods (about 3 hours).
Computational experiments show the efficiency and comparative computational efficiency of the proposed method. When solving problems (4) -(9) by the proposed method, all 6 problems are solved simultaneously and it is also not necessary to solve systems of equations of large dimension, there are no iterative cycles. While solving these problems by continuous optimization methods, each problem is solved separately.

Conclusions
It is shown that the problem of dimension is an obstacle to the creation of methods suitable for optimizing the modes of heat supply systems of large dimension, but at the moment there are no correct methods for solving it.
To solve the dimensionality problem, it is proposed to use a hierarchical approach to optimizing the HSS modes.
The mathematical formulation of the problem of aggregation of distribution heating networks, arising at the stage of planning modes, following from the problem of optimization of the HM HSS and the hierarchical approach to the problem of optimizing the HM HSS.
To solve the problem of DHN aggregation, an original modification of the loop reducing dynamic programming method is proposed. The advantages of this method are: the absence of iterative processes, the absence of the need to solve systems of equations of large dimension, the possibility of solving the problem of DHN aggregation "in one pass".
The proposed method is implemented as a research program. The results of its application for numerical calculations are presented, illustrating its performance and better computational efficiency in comparison with the reference methods.
The study was carried out within the framework of project III. 17