Optimization of hydraulic conditions of radial district heating systems with pumping stations Текст научной статьи по специальности «Математика»

Optimization Of Hydraulic Conditions Of Radial District Heating Systems With Pumping Stations

Alexander Lutsenko, Nikolay Novitsky Melentiev Energy Systems Institute of Siberian Branch of Russian Academy of Sciences, Irkutsk, Russian Federation

Abstract - The paper focuses on the problem of multi-criteria optimization of hydraulic conditions of district heating systems on the basis of the following objective functions: minimum costs of maintaining the operating conditions, the minimum time for adjustment measures, a decrease in the overall pressure level in a network to minimize unproductive expenses and risks of emergency situations. Specific features of the considered problem are a radial configuration of networks, availability of pumping stations and one heat source. The operating conditions optimization problem is mathematically formalized using controls such as a change in the number of operating pumps, impeller speed, impeller diameter at pumping stations and flow throttling. A modification of the problem-solving method is based on a combination of the methods of dynamic programming and reduction of series and parallel connections of the calculation scheme branches. The final algorithm comprises procedures for creation, stepwise increase, rejection and aggregation of possible phase variable trajectories in the first (forward pass) stage that terminates with the determination of the optimal solution for the single branch with subsequent restoration of components of this solution in the second (backward pass) stage of the algorithm for the initial scheme. The case study illustrating the application of the proposed method to a real district heating system with four pumping stations demonstrates the optimality of the obtained solution by several criteria concurrently and high computational efficiency compared to the methods for a general case of multi-loop networks with unknown flow distribution.

Index terms: multi-criteria optimization, ndistrict heating system.

* Corresponding author. E-mail: luc_alex@mail.ru

http://dx.doi.org/10.25729/esr.2019.01.0005 Received October 05, 2018. Revised April 24, 2019. Accepted May 24, 2019. Available online June 25, 2019.

This is an open access article under a Creative Commons Attribution-NonCommercial 4.0 International License.

© 2019 ESI SB RAS and authors. All rights reserved.

I. Introduction

The problem of an increase in the energy efficiency of district heating systems (DHSs) is topical in Russia and abroad [1]-[4]. The natural processes of equipment ageing and wear, the change in load level and structure, organizational division of technologically related DHSs, their untimely reconstruction and adjustment, and other reasons, all lead to off-design conditions of DHSs, high losses of heat carrier and thermal energy, increase in energy consumption for pumping, emergency rates and violations of technological requirements and uninterruptable heat supply to consumers. The optimal operating conditions of DHSs can contribute to significant energy saving.

The problems of operation optimization arise in different stages of decision making on the DHS control: 1) during reconstruction - to assess the effects of old equipment replacement and new equipment installation (pumps with impeller speed control, automatic regulators and others); 2) during preparation for the next heating season in the process of DHS operation (operation planning) - to adjust and regulate heat networks in view of the changes in the DHS structure and parameters; 3) during real-time control - to plan the loading of basic equipment for the coming day. This paper is devoted to optimization problems of steady-state hydraulic conditions of DHSs that arise in the first and second stages.

In practice, the problems of DHS operating condition planning are solved by the multivariate flow distribution calculations to analyze the consequences of possible measures [5]. In this case, such measures are chosen by an expert performing calculations and the solution quality depends on both the expert's experience and skills and the DHS scale and complexity. As a result, such an approach does not guarantee optimal operating conditions and sometimes even feasible operating conditions. Automation of these problems-solving process is difficult because of some factors, such as their high dimensionality (reaching many hundreds of thousands of variables), nonlinearity, discreteness of part of variables, availability of several optimality criteria, etc. For these reasons, there are currently no methods and software for DHS operation optimization that would be appropriate for wide application. Therefore, the development and application of independent methods,

algorithms, and programs are relevant for the calculation of DHS optimal conditions.

Recently, the optimization problems of DHS operating conditions have been much discussed in scientific publications, which indicates the increasing relevance and at the same time, complexity and versatility of such problems. A large number of studies ([6], [7] and others) are limited to the local small DHSs. The techniques proposed for large DHSs are 1) the aggregation of DHS schemes ([8], [9] and others), which do not allow considering the whole set of initial conditions and constraints; 2) the approximation of relationships between the objective function value and the state variables (chosen as the basis) [10], [11] in case of difficulties in adequately considering discreteness of the optimized equipment configuration; 3) the application of the off-the-shelf software to solve optimization and simulation problems by the general-purpose mathematical methods, which do not take into consideration the DHS features, and require considerable computational effort. In [12], [13], for example, the model of a physical network is constructed in the Simulink/Matlab environment and its state is calculated by the off-the-shelf CPLEX solver, while its optimization is performed by the ReMIND software; 4) the application of the semi-heuristic, but easily-implementable genetic and evolutionary algorithms, which, however, involve even more computational effort. In [14], for example, the pumping expenses are minimized by the nested iteration cycle. In the internal cycle, the feasible condition is iteratively calculated by the SIMPLE-algorithm. In the external cycle, the parameters that influence the condition are changed using the genetic algorithm. Optimization of the DHS conditions often involves partial criteria. For example, in [15] the total fuel consumption is minimized at a specified heat carrier flow rate in the pumps. In [10], [11] the optimal control problem of pumping stations (PSs) is solved using the criterion of minimum consumed electricity under specified heat carrier temperatures. At the same time, there are practically no works devoted to the multi-criteria statements of the optimization problems or the multi-criteria problem is solved by the formal weighting of diverse criteria [16].

In addition, the dynamic programming (DP) method is applied effectively enough and widely to solve the optimization problems of expansion and reconstruction of the tree-structured DHSs (when the flow distribution is known). This method is potentially applicable to the optimization of operating conditions ([17], [18] and others) as well. For a general case of the multi-loop networks with the unknown flow distribution, S.V. Sumarokov suggested a method based on the iteration process with the alternate application of the methods of flow distribution calculation and DP [19], [20]. A modification of this method [21], which also involves the iterative process of solution search but has no rigorous substantiation was proposed for multiloop pipeline systems with the specified flow distribution.

There are also modifications of the DP algorithms

for optimization of radial DHSs. These modifications are based on separate optimization of the supply and return pipelines of a heat network [18]. The radial DHSs have a tree topology in a single-line representation but their calculation schemes are multi-loop ones because each consumer forms an independent loop of heat carrier circulation. Applicability of such decomposition to the optimization of operating conditions was studied in [22]. The study indicated that solving the problem "by parts" did not yield either an optimal or even a feasible solution for the whole DHS.

This paper presents the results of the evolution of the DP method modification suggested by the authors earlier for optimization of hydraulic conditions of distribution heat networks [22] in order to allow for active components (pumping stations) and several optimality criteria. The proposed DP modification is applied to hydraulically related DHSs of a radial structure with one heat source (HS) which are studied as a whole (without decomposition). The temperature curves at a heat source are assumed to be specified, the heat losses are eliminated, and their residual value may be neglected. In this case, the requirements of heat energy supply to consumers are reduced to the necessity of maintaining an appropriate heat carrier flow rate, and the problem is reduced to hydraulic conditions optimization. Fixing the flow rates in the radial DHSs specifies the unique flow distribution in the network, and therefore the subject of DP method study and evolution is the hydraulic condition optimization in the network with the specified flow distribution and the multi-loop topology of a special type. Provided there are pumping stations in DHS, they are taken into account through consideration of the case typical of a DHS when the similar centrifugal pumps are installed at the pumping station and connected in parallel.

The proposed DP method with loop reduction [22], [23] has some advantages: 1) linear increase in the computational effort depending on the problem dimension, which allows optimization of large DHSs; 2) verifiability of the existence of at least one feasible solution; 3) possible optimization on the basis of several criteria simultaneously; 4) guaranteed optimality of the obtained solution; 5) high computational efficiency (without the iteration process) compared to the methods for a general case of multi-loop networks with the unknown flow distribution [24]; 6) potential adaptability to consideration of pumping stations while retaining all these properties.

II. Problem statement

The substantive objective of optimizing the hydraulic conditions of DHS is to find controls that ensure the implementation of a condition that satisfies the feasibility and optimality requirements in accordance with a given system of criteria.

The feasibility requirements are reduced to the necessity of meeting consumer loads and observing conditions of

equipment operation. The energy saving requirements may be reduced to the uniform economic criterion - the variable component of operating costs. The other (technological) criteria are intended to minimize labor input of adjustment measures, decrease heat carrier leakages and risks of emergency situations. Optimization on the basis oftechnical criteria may be reduced to the minimization of pressure control points and minimization of the overall pressure level in the network. The considered controls are aimed at a possible change in the performance of DHS components (pumping stations, pipeline sections, and consumers) and include a change in the number of operating pumps, the impeller speed, the impeller diameter, and flow throttling. For the pipelines and the sections modeling consumers, the controls are performed mainly by throttling (by the regulator, the balancing valve or other devices).

III. Models of the controlled DHS components

Denote sets of the indices of branches of the DHS calculation scheme that model pipelines, pumping stations, and consumers by IPL, IPS, IC, respectively. Then IPL u IPS u IC = I is the set of indices of all branches,

|/| = n , n is the number of the scheme branches.

A generalized hydraulic characteristic of the branch can be represented by [24]

h. (x., z. ,y. ,k. ) = zs.x.lx.lIk2-v2H., i e I (1)

iii|i| i ' i i7

Here: i is the branch index; hi is the pressure drop; xi is the heat carrier flow rate; si is the hydraulic resistance of the pumping station; Hi is the pressure rise caused by pumps at the pumping station; yi is the relative impeller speed (or its diameter); ki is the number of operating (from installed) pumps at the pumping station; zt is the hydraulic resistance index because of throttling (z. > 1). Thus, z, yi are the continuous controls, and k are the discrete controls.

In (1) any type of unavailable or impermissible control can be taken into consideration by assigning the constant

to its value. For example, for i e IPL u IC : Hi = 0 , ki = 1.

For the illegal throttling zt = 1, i e I . If the speed control is impossible, then yi = 1, i e Ips. The variant ki = 0, i e Ips in (1) is modeled by the values ki = 1, Hi = 0, and si is the bypass line resistance with possible throttle installation when zt is variable.

Model of controlled hydraulic condition. Based on the indicated features of the radial network topology (the tree topology in a single line representation, the multiloop calculation scheme, symmetric connection schemes of components for heat carrier supply and return which are connected through the source and consumers) (Fig. 1), and the requirement to meet the specified flow rates at consumers, we have fixed distribution of flows in the network (x, ie I). Then expression (1) can be substituted by

(z. ,y. ,k. ) = z.s.x\x\ Ik2 -y2 H., i e I, and the model of

* 1 ' 1 i' iii|i| i ' i i7

the controlled hydraulic condition of DHS which is applied in [24] for the general case of unknown flow distribution will have the form

f

A P - y y - œ(z,y,k)

= 0

(2)

where: A is the rn^n-dimensional incidence matrix of the DHS calculation scheme;

m is the number of nodes in the calculation scheme; Q, P are the m-dimensional vectors of nodal flow rates and pressures; y is the n-dimensional vector of pressure drops in branches; m(z, y, k) is the n-dimensional vector function with the components

&r(z,,Yi), i = 1,n .

1

10

I

—7

O^4

z1

s I As a

s

(b)

0

a

3—

m

Fig. 1. Example of a single-line technological (a) and double-line calculation (b) DHS schemes and a block diagram of pumping station (c). 1 - heat source, 2 - heat chamber, 3 - consumer, 4 - pumping station, 5 - pipeline, 6 - control point z,, 7 - heaters, 8 -pump, 9 - bypass line, 10 - number of node j (the figures in circles), 11 -number of branch i (the figures in squares).

Let us call part of the hydraulic condition parameters that depend on the environment boundary conditions. When the heat source is a boiler plant, it can be treated as an ordinary pumping station with the given pressure at the makeup node (at the inlet). Optimization of the cogeneration plant operation is an independent subproblem [26], and consideration of its operation at known connected load is reduced to assigning two pressures (at the inlet and outlet).

The feasibility requirements are: p < p < p , j = 1,m,

y < y < y , i = 1, n, where the bottom and top lines

indicate the given boundaries for the feasible value of the same parameter. The requirement of the effective range of feasible capacities should also be observed for pumping

stations. Let X' X be the boundaries of the effective range of flow rates for one typical pump on the z-th branch with a pumping station, then we have the inequality

KYiXi < x. <Kiyixi or xj% < K.ji < xjXi i e Ips.

Besides, it is necessary to comply with the feasibility requirements of the controls: 1 < z. < z., y. < y. < y.,

ki e {0,1,2,...Kz}, z e I. Introduce the vector of the Boolean variables 8, whose components 8ze{0,1} indicate the absence or presence of a throttle on the z-th branch. The inequality 1 < z. < z., in this case, may be expressed as

1 < z. < 1+(z. -1) 8. On the whole, the system of inequality constraints has the form

P < P < P , y. < y. < y.,

—j j j > Si Si?

1 < Z,. < 1+(z -1)5,., Y <yj < Y, (3)

K e(0,l,2...K,}, xj Xi ± KYi < xj % , i e I

Objective functions. The variable component of costs of maintaining the operating conditions consists of electricity costs (to pump heat carrier) and fuel costs (to heat water). In the considered case of DHS with one heat source, its heat load is fixed (equal to a total load of connected consumers), the fuel costs do not depend on the heat network operation and can be optimized within the local problem at the level of the heat source scheme. The electric power consumed by the pumping station is set by the function

N (YK) = + + ^YSt'K,, * * trs [24].

Here fi0, are the coefficients of the square polynomial that approximates the power consumption characteristic by one pump. The economic objective function has the form Fc (y,k) = X cEPNi (y,,Ki), where cEP is the electricity

ielps

rate for the i-th pumping station, y,k are the real and integer vectors with the elements y, k,, z e Ips. Correspondingly, the "technological" objective functions [25] are the number of control points - Fz (5) = ^8. ; the average pressure in the

network - Fp (P) = ^ P. / m, where J is the set of indices of all nodes. JeJ

IV Mathematical problem statement.

The known are: the calculation scheme topology (matrix A); the fixed pressure at least at one node (P*); the coefficients of hydraulic and capacity characteristics (Po , Pi,, P2 ,, , e Ips, s , H , z e I); the acceptable limits of change in continuous parameters of conditions and controls; the sets of possible values of the integer parameters ({0,1,2,... K}, z e Ips ); and the electricity cost ( cf,ie Ips ).

By solving the following optimization problem with F=(Fc,Fz,Fp)T min F subject to (2), (3), (4) it is necessary to determine the throttling sites (8) and value (z), the number of operating pumps (k) and the values of their speed (y), as well as the operating parameters (P and y).

The great number of criteria can be taken into consideration based on the lexicographic ordering of objective functions by the importance that corresponds to the sequence of their listing in (4). The solution to (4) in this case conforms to the solution to the minimization problem of Fp(P) subject to (2), (3), FC(y,k)=F,* FZ(8)=FZ* where FZ* - is the optimum objective function value in the minimization problem of FZ(8) subject to (2), (3), FC ( y,k)=Fc*, and Fc*- is the optimum objective function

value in the minimization problem of FC(y,x) subject to (2), (3).

V. Problem-solving method

Basic concepts. In terms of the DP theory, the pressure (Pj) can be determined as a phase variable concentrated at the calculation scheme nodes, and the pressure drop (yz) - as the main parameter of the elementary trajectory section of this variable along one branch, the remaining unknowns (z,y,x) - as controls of this trajectory. DP involves 1) the possibility of representing the process of solution search as a multistep process; 2) the independence of decision making at each step of the stepwise increase in the phase variable trajectory of its prehistory; 3) additivity of the objective function that serves to verify its value for any section of this trajectory. A multistep process of the stepwise increase in the phase variable trajectories can be implemented by bypassing the scheme branches through the nodes connecting them. The only simple path connects any two nodes in the tree networks. In the considered case with available loops, the things are quite different, and the main problem is to reject trajectories that do not satisfy the second Kirchhoff law. The second condition is observed owing to the fixed flow distribution since the phase variable trajectory at each step depends only on the parameters of the function mi(z.,y.,k.) at this step. All of the above objective functions are additive.

The main idea of the proposed modification of the DPLR (DP with loops reduction) method consists in a combination of the methods of DP and reduction of the series-parallel connections of the calculation scheme branches. The main DP procedures include creation,

stepwise increase, and rejection of possible phase variable trajectories in the forward pass that terminates when the optimal solution is obtained, with subsequent restoration of the solution components in the backward pass of the algorithm. The reduction techniques are based on two simple relations ^e ~ X y and y„ = y , i e R for the

f s E S,' 2

pressure drop in the equivalent branch (ye) that substitutes the set R of branches connected in series or the set R2 of branches connected in parallel. The last relation follows from the second Kirchhoff law for the simplest loop of two branches.

Trajectories of the phase variable are created in the space of its discrete values. The feasible pressure range

at each node is divided into W. = (P - P )/ e„ cells,

j j j p '

where is the specified cell size. Therefore, Wj possible discrete pressure values P (w ), w = 1,W correspond to each node. Denote the starting and end nodes of the i-th branch by f e J and l. e J , respectively. By varying the feasible controls on this branch for each Pfi (wf .), wf. = 1, Wfi at the node f, we generate various trajectory sections such that yk = (zk ,yk. ,Kk), y < yk < yi and

№,) = Pf,(wfl)-y eiZ.Mi where Pn(wh) is rounded to the value corresponding to the nearest cell center at node l ,, and the unique combination of wH, wf. corresponds to each k. Thus, the cells (possibly, not all) are filled at node l , and on each branch there emerges a set of feasible local trajectories of the phase variable with pressure drops, which form the set Y = U {yk} , where k is the index of the element of this set. Each trajectory section

on the -th branch can be associated with the increase Fk

i

in the value of the vector objective function F with the components: (FC )k = of N. {jkt ,Kk); (F )k = S, where Sk = 0 at zk = 1 and Sk = 1 at zk > 1; (F )k = Pk, where

i i i i V P'l ll

Pk is the pressure at the end node for yk. In the case of a pumping station, there can be a non-unique combination of controls (and correspondingly, of the values (F )k and (Fz )*) that implement yk. At the same time, for the i e IPL u IC control z e[1, z ] uniquely corresponds to the specific y . Therefore, yki e IPS should be generated using the local procedure for selection of the optimal combination of controls.

vI. Local optimization of a pumping station Local optimization of a pumping station is aimed at determining the controls zk ,jk ,Kk optimal by the criterion F , which implement the given , ' e 1PS .

For pumping stations, without speed control, the problem is reduced to the enumeration of the options k = 1,2,... and for each option z e [1, z ] is determined and FK is

calculated (subject to the condition x./- xi/X)

based on the relation yk = (z ,yk,k ) (where yk = 1).

For the case of speed control (presuming that it consumes less power than throttling), the problem is also solved by enumeration of the options k = 1,2,...,K , and for each: 1) the solution existence with respect to y i which is reduced to verifying the condition < fK. , where

YK = max,xt/(xk, )], Tt = min\_fi,x,/(xk)] , is tested; 2) the solution y* to the equation yk = & (z',y ,k ) where z* = 1, is found; 3) if ^ Y* ^ YK , the solution z*,j* is found; if y* < yK , then assuming that 7* = yk , z*

11 1 — 1 r i t_ j i

is defined as the solution to the equationyk = mi(z. ,y',k. )., Otherwise, the solution z*,y* does not exist, because the increase in throttling causes the increase in y., and that of speed control - to its decrease; 4) if the solution is found, FK is calculated.

Whether or not there is an speed control: 1) the option k. = 0, which is reduced to identifying z. e [1, z. ] to implement yk, is additionally tested (just as for i e IPL u IC . If the feasible solution z* exists, then FK

is determined; 2) if there are no feasible solutions among the considered variants, then the k-th trajectory section is rejected; if there are several of them, the solution zk ,Yk ,Kk that corresponds to Fk = minFK is chosen.

vil Extension and aggregation of trajectories

Consider two series-connected branches that have

common node l.. = f.„ and no other branches are

i\ j ii

incident to this node ln = f.2. Let the set Y.1 = U {>"*} be formed for branch i1, therefore the cells for the pressures Pfi2(wfi2) are filled. Then the set Yi2 = U{yk2} is generated from these cells, and, as a result, the trajectories

pfn(wfIi) - pm(wM) - pn 2(wu 2)=$+yk:, where

P,n(wM) = Pfl2(Wfi2), are built.

For the equivalent branch i 3(substituting i 1 and i 2) such

a trajectory is aggregated by the rules: y.3 = ya + y

(Fc )k3 = (Fc + (Fc )k?2 , (Fz)» = (Fz)» + (Fz 12 ,

(Fp )k3 = (Fp)» + . If in the process of such aggregation there occur the trajectories with the coinciding end pressures (wfi3 = w = wu 2 = wu 3), the resulting trajectory (k3) is

given the best value of increase in the vector criterion F3. Eventually, the set Y.3 = U{yk3} is formed.

Let there be two branches (i1, i2) connected in parallel,

such that f = f , 1 = 1 , and for one of them (for

Ji 1 Ji2 ¡1 12' V

example, i1) the set Yn is generated. The coincidence of the end pressures (w n = w ,2 and wln = wH2) for each pair

y", yi2 means that they satisfy the second Kirchhoff law. Therefore, for branch z2, it remains to verify feasibility of such trajectories that y^ = y" e Yi1. The feasible pairs of trajectories that correspond to yyk^ are aggregated into the trajectory of the equivalent branch z'3 (substituting z1 and z'2)

by the rules: = yi = (Fc£ = (Fc+ (Fc fl

(F £ = (F tl + (F (FP )"a = (F tl = (F

The remaining trajectories are not considered. The set Y.3 = U{yk3} is formed as a result of such processing

VIII. The computation scheme of DPLR, in general terms, includes the following basic stages.

1. Select the next fragment of the series-connected branches. If there are no such fragments, then go to stage 3. In the initial scheme of the radial DHS such fragments are always available. At least, any dead end with a consumer forms a series connection of three branches in the calculation scheme (for the supply line, the consumer, and the return line, see Fig.1).

2. Replace this fragment with one branch using the considered procedures of generation, extension, rejection and reduction of the phase variable trajectories for series connections of the branches. If in the generation process Y. = 0 even for one initial branch, then go to stage 6. If the current DHS scheme

Fig. 2. An example of reducing a calculated scheme of a heatingsystemtoone branch byreducing seriesand parallel connections of branches. In the calculated scheme: bold lines -seriesconnections of branches, dottedlines - parallel connections.

contains only one branch, then go to stage 5.

3. Select the next fragment containing the branches connected in parallel. As a rule, at least one of such branches equivalents series connections processed in the previous stage. Therefore, replace this fragment with one branch using the considered aggregation techniques of the trajectories of parallel connections. If in the aggregation process Y = 0 at least for one branch, then go to stage 6.

4. If the current DHS scheme contains more than one branch, then go to stage1.

5. The trajectory of the optimal value of Fk for the only remaining branch corresponds to the aggregated problem solution. The backward pass of the algorithm is used to restore the optimal phase variable trajectory and the values of all controls in the initial scheme. This suggests that the correspondence between each equivalent branch and the substituted fragment and the correspondence between each section of the equivalent branch trajectory and the initial sections it aggregates are stored in the previous stages of the algorithm. End.

6. There is no feasible problem solution.

This algorithm is intended for the radial DHS with at least one heat source and one consumer that are connected by the branches for the supply and return lines. The reduction process of the radial DHS scheme to one branch is illustrated in Fig. 2. The Figure shows that there are series connections in the initial scheme, and their equivalenting (reduction) results in the formation of parallel connections, which in turn leads to the formation of new series connections and so on until the calculation scheme is reduced to one branch connecting the nodes at the heat source inlet and outlet.

Ix. Numerical example

The proposed method was tested by a set of computational experiments. Their results are illustrated by the example of DHS in the town of Baikalsk, which is presented in an aggregated form in Fig. 3. PS-1 is installed on the supply pipeline, and the remaining PSs - on the return pipeline. Table 1 describes the characteristics of the installed pumps. All pumps are supplied with electricity at the same cost, therefore FC (j,k) = ^ N. (y. ,k. ).

PS

Reliability and computation efficiency of the proposed DPLR method were tested in comparison with the continuous branch and bound method (CBBM) [24] that was developed for a general case of the discrete-continuous optimization problem of hydraulic conditions in multiloop DHSs with several heat sources, pumping stations, unfixed flow distribution but with one criterion (FC). The calculation results by CBBM coincided with DPLR, however, in the last case, the computing time was less by anorderof magmtudo.

Step Transformation of the calculated (two-line) scheme Transformation technological (single-line) scheme

1 A h Sï v

2 j t V 1 1/ ' \ %

3 l\_

4 * '*» t' /' D V V

5 3

6 1 N

Fig. 3. Aggregated scheme of DHS in Baikalsk Pumping stations on the return pipeline are highlighted in grey.

Besides, as was shown above, the computing time by DPLR did not depend on the number of criteria, whose lexicographic ordering is performed by elementary operations at each step of the DP process. Therefore, the multi-criteria problem is solved in a single algorithm run. Such problems are solved by the most common approach that involves sequential optimization based on multi-stage solving the optimization problems for each criterion by its importance with intermediate addition of constraint on the next criterion value calculated in the previous stage. Therefore, considering the great number of criteria by the standard methods in this case will additionally increase the computing time at least threefold, and the mixed nature of the problem (nonlinear, integer, Boolean programming) will increase it by several times more.

The flexibility of the proposed DPLR method (consideration of several criteria, availability of pumping stations and different control techniques) was tested for DHS in Fig. 3 for different combinations of calculation data: 1) different levels of consumer loads including "heavy" load (corresponding to the peak winter load of DHS) and "light" load (obtained from the previous one by the twofold load decrease); 2) impossibility (y. = const = 1, i e IPS) or possibility (y = var, i e I ) of pump speed control at the pumping stations.

Table 2 presents the results of optimization calculations in comparison with some feasible condition ("heavy feasible") that meets all constraints. The obtained solutions are seen to be more efficient for all criteria compared to this condition. Besides, the speed control considerably decreases the power consumed by pumping stations and the number of sites for throttling. The twofold load reduction in DHS ("light" load) leads to a more than threefold decrease in the power consumption, in particular owing to a decrease in the number of operating pumps and their complete disconnection at PS-3. A decline in the overall pressure level in DHS is notably lower than for an arbitrary feasible condition and it can be additionally decreased by the lower heat load.

X. conclusions

1. For the first time, the problem of operating condition planning for the radial DHSs with pumping stations was formulated and studied as a mixed multi-criteria optimization problem with real, integer and Boolean variables. At present, there are no effective methods appropriate for extensive practical application.

2. The paper proposes a generalized option of the problemsolving method which is based on a combination of the methods of dynamic programming, reduction of seriesparallel connections of branches in the calculation scheme and local optimization of pumping station operation.

3. The presented numerical example illustrates the computational efficiency of the proposed method and its versatility, which makes it possible to allow for one or several criteria, different types of controls and technological limitations, arbitrary number and allocation of pumping stations, establishment of solution existence, applicability to multilevel coordination of optimal conditions for heat networks and heat sources.

References

Table 1. Characteristics of pumps at pumping stations.

PS Number of pumps (K) Pump brand H s ßo ß ß

1 3 SE-800-100-11 120 0.0000375 119.4 0.238 -0.000091

2 5 1X 200-150-500 80 0.00004 40 0.16 0

3 3 TsN-400-105 120 0.000125 60 0.24 0

4 2 D320-50 60 0.0001 30 0.09 0

Table 2. Results of optimization based on the number of operating pumps at pumping station.

Load yi > ie ips PS-1 PS-2 PS-3 PS-4 FC, kW FP m ofthrwatercolumn

pcs.

Heavy feasible 1 3 4 3 1 1684 20 62

Heavy 1 3 4 2 0 1445 10 58

Heavy Var 3 4 2 0 1396 7 58

Light 1 1 2 0 1 426 10 56

Light Var 1 2 0 1 409 7 56

[1] I.A. Bashnakov, "Analysis of primary trends in heat supply system development in Russia," Novosti Teplosnabzheniya (News of Heat Supply), No. 2-4, 2008, (in Russian.)

[2] V. Roshchanka, M. Evans, "Playing Hot and Cold: How Can Russian Heat Policy Find Its Way Toward Energy Efficiency?", National Technical Information Service, USA, p. 30, 2010.

[3] Henrik Lund, Sven Werner, Robin Wiltshire, Svend Svendsen, Jan Eric Thorsen, Frede Hvelplund, Brian Vad Mathiesen, "4th Generation District Heating (4GDH) Integrating smart thermal grids into future sustainable energy systems", Energy, № 68, p. 1-11, 2014.

[4] "Energy Policies of IEA Countries," Canada, 2009, Review, OECD/IEA, 2010.

[5] V.V. Tokarev, Z.I. Shalaginova, "Technique of Multilevel Adjustment Calculation of the Heatthydraulic Mode of the Major Heat Supply Systems with the Intermediate Control Stages", Thermal Engineering, vol. 63, № 1, p. 68-77, 2016, DOI 10.1134/S0040601516010110.

[6] D. Buoro, P. Pinamonti, M. Reini, "Optimization of a Distributed Cogeneration System with solar district heating," Applied Energy, vol. 124, p.p. 298-308, 2014, DOI: 10.1016/j.apenergy.2014.02.062

[7] M. Larsson, J. Dahl, "Reduction of the specific energy use in an integrated steel plant - the effect of an optimization model". ISIJ Int; 43:1664-73, 2003.

[8] H. Larsen, "A comparison of aggregated models for simulation and operational optimization of district heating networks", Energy Convers Manage, 45: 1119-39, 2004.

[9] H. Larsen. "Aggregated dynamic simulation model of district heating networks," Energy Convers Manage, 43: 995-1019, 2002.

[10] S. Cosentino, E. Guelpa, R. Melli, A. Sciacovelli, E. Sciubba, C. Toro, V Verda, "Optimal operation and sensitivity analysis of a large district heating network through POD modeling," Proceedings of the ASME 2014 International Mechanical Engineering Congress and Exposition, Montreal, Quebec, Canada, 2014.

[11] A. Sciacovelli, E. Guelpa, V. Verda. "Pumping Cost Minimization in an Existing District Heating Network," ASME 2013 International Mechanical Engineering Congress and Exposition, vol. 6A: Energy, p. V06AT07A066, Nov. 2013, DOI:10.1115/ IMECE2013-65169

[12] Elisa Guelpa, Claudia Toro, Adriano Sciacovelli, Roberto Melli, Enrico Sciubba, Vittorio Verda, "Optimal operation of large district heating networks through fast fluid-dynamic simulation," Energy, vol. 102, p.p. 586-595, 2016, ISSN 0360-5442, DOI

10.1016/j.energy.2016.02.058.

[13] Mattias Vesterlund, Jan Dahl. "A method for the simulation and optimization of district heating systems with meshed networks", Energy Conversion and Management, 89, p. 555-567, 2015.

[14] Adriano Sciacovelli, Elisa Guelpa, Vittorio Verda, "Pumping cost minimization in an existing district heating network", Proceedings of the ASME 2013 International Mechanical Engineering Congress and Exposition IMECE2013, November 15-21, 2013, San Diego, California, USA IMECE, 65169, 2013.

[15] X.S. Jiang, Z.X. Jing, Y.Z. Li, Q.H. Wu, W.H. Tang, "Modelling and operation optimization of an integrated energy based direct district water-heating system", Energy, vol. 64, Jan. 2013, DOI. 10.1016/j. energy.2013.10.067.

[16] Hai Wang, Haiying Wang*, Zhou Haijian, Tong Zhu, "Optimization modeling for the smart operation of multi-source district heating with distributed variablespeed pumps", Energy, vol. 138, p. 1247-1262, 2017.

[17] A.P. Merenkov, V.Ya. Khasilev, "Theory of hydraulic circuits", Moscow, Nauka (Science) Publishers, 1985, (in Russian.)

[18] E.V. Sennova, V.G. Sidler, "Mathematical modeling and optimizing developing heat-supplying systems", Novosibirsk, Nauka (Science) Publishers, 1987, (in Russian.)

[19] S.V. Sumarokov, "Method to solve an extremal network problem", Ekonomika i matematicheskiye metody, (Economy and mathematical methods), vol. 12, No. 5, p.p. 1016-1018, 1976, (in Russian.)

[20] A.P. Merenkov, T.B. Oshchepkova, S.V Sumarokov, V.Ya. Khasilev, A.V. Khramov, "Optimal synthesis of multi-loop systems with loaded reservation", In: Energy Systems - trends of development and methods of control. Irkutsk: SEI SB AS USSR, vol. 1, p.p. 180-192, 1980, (in Russian.)

[21] A.P. Merenkov, E.V Sennova, S.V. Sumarokov, et al. "Mathematical modeling and optimizing systems of heat-, water- oil-, and gas supply", Novosibirsk, Nauka (Science) Div., Siberian Publishing Firm, 1992, in Russian.

[22] A.V. Lutsenko, N.N. Novitsky, "Modified method of dynamic programming for optimization of hydraulic modes of distribution heating networks",

Computational Technologies, vol. 23, No. 6, p.p. 48-64, 2018, (in Russian) DOI: 10.25743/ ICT.2018.23.6.006.

[23] A.V. Lutsenko ""Applying the dynamic programming method at multilevel optimization of heat-supplying system hydraulic modes," Sistemnyye issledovaniya v energetike (System research in power engineering). Transactions of young scientists at MESI SB RAS, vol.

47, p.p. 32-31, 2017.

[24] N.N. Novitsky, A.V. Lutsenko, "Discrete-continuous optimization of heat network operating conditions in parallel operation of similar pumps at pumping stations," J. Global Optim, vol. 66, No.1, p. 83-94, 2016, DOI. 10.1007/s10898-016-0403-y.

[25] A.V Lutsenko, N.N. Novitsky, "Investigating problems and methods for multi-criterion optimization of supply heating network hydraulic modes," Scz. Bull. NSTU, No. 3, vol. 64, p.p. 131-146, Jul. 2016, DOI. 10.17212/1814-1196-2016-3-131-145.

[26] "Effective methods for optimization of schemes and parameters of heat power facilities: development and application" / Ed. by A.M. Kler. Novosibirsk:

Academzc PubHshzng House "Geo", p.210., 2018, ISBN (print): 978-5-6041446-0-2.

Lutsenko Alexander Viktorovich,

senior engineer of the Institute of Energy Systems of the Siberian Branch of the Russian Academy of Sciences (Irkutsk). His research interests are optimization methods and their application to pipeline systems. Author of more than 15 published scientific papers on this subject.

Novitsky Nikolay Nikolaevich, Doctor of Technical Sciences, chief researcher at the Institute of Energy Systems of the Siberian Branch of the Russian Academy of Sciences (Irkutsk). Area of scientific interests - the theory and methods of analysis, management, optimization and identification of pipeline systems; mathematical and computer modeling. Author of more than 200 published scientific papers on this subject. Chairman of the All-Russian permanent scientific seminar with international participation "Mathematical models and methods for analysis and optimal synthesis of pipeline and hydraulic systems." In recent years, he has been actively engaged in theoretical and applied problems of the intellectualization of pipeline systems.