A Nonlinear Analytical Algorithm for Predicting the Probabilistic Mass Flow of a Radial District Heating Network

This paper develops a nonlinear analytical algorithm for predicting the probabilistic mass flow of radial district heating networks based on the principle of heat transfer and basic pipe network theory. The use of a nonlinear mass flow model provides more accurate probabilistic operation information for district heating networks with stochastic heat demands than existing probabilistic power flow analytical algorithms based on a linear mass flow model. Moreover, the computation is efficient because our approach does not require repeated nonlinear mass flow calculations. Test results on a 23-node district heating network case indicate that the proposed approach provides an accurate and efficient estimation of probabilistic operation conditions.


Introduction
The development of a low-carbon sustainable energy system has generated increasing interest since energy and environmental issues have become more prominent globally.The present direction of this development has focused on integrated energy systems (IESs) that integrate various energy-related tasks, such as cooling and heating, and various forms of energy, such as electricity and natural gas, to provide a comprehensive utilization and management of energy [1].However, the increasing integration of energy systems, such as combined heat and power (CHP), gas turbines, and other energy conversion facilities, has greatly increased their interdependence [2,3].This interdependence necessitates increasingly sophisticated planning and operation of energy systems.
Presently, the planning and operation of IESs is generally based on the steady-state modeling and analysis of IESs.For example, Liu et al. [4,5] established a steady-state model between electricity and heating networks and proposed an effective mass flow calculation method.Similarly, steady-state energy-flow analysis between electricity and natural gas networks has been conducted [6,7].Moreover, steady-state energy flow analysis has been conducted with integrated electricity-gas-heat systems [8][9][10].However, steady-state analyses are essentially deterministic and cannot effectively address the many uncertainties arising in IESs, such as random fluctuations in cooling, heating, electricity, and natural gas loads; intermittent energy output fluctuations; generator failures; electric transmission line and gas pipeline failures; and market uncertainties.Moreover, different energy networks have different mutual influences affecting their interdependencies.These studies highlight the necessity of investigating the planning and operation problem of gas/heat systems that are interdependent with power systems.
Recently, numerous studies have sought to develop analysis methods capable of addressing the influence of uncertainties on power networks [11][12][13].For example, Chen et al. [14] considered the uncertainties of electricity, gas, and heat loads and wind farm outputs within the framework of steady-state energy flows and employed Monte Carlo simulations to solve the probabilistic energy flow of IESs.The effects of various uncertainties on IES reliability have also been studied [15].A probabilistic steady-state analysis of integrated electricity, gas, and heating networks has been proposed based on Latin hypercube sampling and the Nataf transformation [16].A stochastic scheduling model is proposed for the interconnected EHs considering integrated demand response (DR) and wind variation [17].In addition, ensemble prediction systems (EPSs) have laid a foundation for the quantitative analysis and evaluation of the influence of uncertain factors in IESs.In an EPS, the statistical characteristics of random output variables can be obtained according to the statistical characteristics of random input variables by the calculation of probabilistic flows.The methods for solving probabilistic flows include simulation methods [18,19], analytical methods [20,21], and approximation methods [22][23][24].Monte Carlo simulation is the most common simulation method employed to test the accuracy of probabilistic power flows.The middle semi-invariant method is the most widely employed analytical method owing to its high calculation efficiency.Finally, the most representative point estimation method represents a very commonly adopted approximation method because it requires no knowledge regarding the specific functional relationship between input quantity and output quantity.At present, the probabilistic power flow calculation in power systems has been extensively investigated.However, it should be clarified that few studies have investigated probabilistic mass flows in thermal systems.Moreover, usage of nonlinear analytical algorithms for solving probabilistic mass flows has not been discussed.
The present work addresses these deficiencies in past works by developing a nonlinear analytical algorithm for predicting the probabilistic mass flow of a radial district heating network based on the principle of heat transfer and basic pipe network theory.First, the variance of mass flow through a pipe connected with a heat source is obtained according to the power balance equation of a district heating network.Then, the functional relationship between the mass flow variances between pipes in the network is deduced to obtain the variance of mass flow in the entire pipe network.Second, a functional expression of the pipe network node temperature is derived, and the covariance matrix of the mass flow through the pipe network is obtained.Finally, the variance of the node temperature can be obtained.The validity and rationality of the proposed algorithm is verified by application to a 23-node radial district heating network with various pipe lengths under thermal load fluctuations of various magnitudes.The probabilistic results obtained can provide comprehensive information of real-time IES operating conditions, which is valuable for IES planning, operations, and risk assessment.

Steady-State District Heating Network Model
The steady-state model of a district heating network includes a hydraulic model and a thermodynamic model.The hydraulic model includes the joint flow equilibrium equation and the pressure head loss equation, which are, respectively, expressed as follows: where, A is the node-branch incidence matrix.m and mq represents, respectively, vectors of mass flow rate within each pipe and the injected mass flow at the nodes (kg/s).hf represents the vector of head losses (m), and K represents the vector of resistance coefficient of pipes.The thermodynamic model includes the thermal load power equation, the pipe temperature change equation, and the node power conservation equation, which are, respectively, given as follows: where Φ represents the vector of the heat power consumed or supplied (MW).Cp is the specific heat of water, and Cp = 4182 × 10 −3 MJ•kg −1 •°C −1 .Ts represents the vector of supply temperatures at nodes (°C).To represents the vector of the outlet temperature of flow at the outlet of nodes before mixing in the return network (°C).Tstart and Tend represents the temperatures at the start node and end node of the pipe, respectively (°C).h represents the total heat transfer coefficient per unit length (W/(m•k)).L represents the length of the pipe (m).Ta represents the ambient temperature (°C).mout and min are, respectively, the mass flow rate leaving and entering the node (kg/s).Tout represents the mixture temperature at the node (°C), and Tin represents the temperature of mass flow entering the mixing node at the end of the incoming pipe (°C).
Equations ( 1)-( 5) are nonlinear equations, where the coupling relationship between temperature and mass flow is strong, and exponential terms are involved.Therefore, solving these equations for realistic thermal pipe networks directly is quite difficult owing to the high computational complexity involved and the inability for ensuring numerical stability.

Probabilistic Thermal Load Model
In general, thermal loads can be described probabilistically in terms of a normal distribution.Accordingly, the thermal load probability density function (PDF) can be described as: where μ and σ are the respective mean and standard deviation of thermal load .

Approximate Model of Probabilistic Mass Flow in a Radial District Heating Network
An approximate model of probabilistic mass flow in a radial district heating network consisting of a heat source node H, three pipes, and three nodes is shown in Figure 1.Here, the circled values represent pipes, and the arrows represent the direction of mass flow rates.As demonstrated in Appendix A, the heat loss of a pipe can be estimated as follows: ( ) The thermal power balance equation for the heat source node H and a pipe network consisting of N pipes and N nodes is described as: where m1 is the mass flow rate of pipe 1, TH is the temperature of the CHP source, To is the return temperature of the CHP source,  i is the thermal load of node i, and   i is the thermal power loss of pipe i. Reordering Equation (8) yields an expression for m1: , where X and Y are statistically independent, then the sum of X and Y also satisfies a normal distribution, i.e., It can be seen from Equation (7) that Δ is approximately constant, so that its variance is approximately zero.Assuming that the thermal load obeys an independent normal distribution, it is known from Lemma 2 that ( ) also obeys a normal distribution with a standard deviation σ.Therefore, the standard deviation of the mass flow rate of pipe 1 can be obtained by Lemma 1 as follows: ( , ) ~( , , , , ) , where ρ is the correlation coefficient between random variables X and Y, then any non-zero linear combination of X and Y also lies within a normal distribution, i.e., For two-dimensional random variables, independence and irrelevance are equivalent characteristics.
The correlation coefficient between mass flow rates can be investigated according to Lemma 3 based on the schematic presented in Figure 2. Here, the correlation coefficient between mass flow rates mi and mj is assumed to be ρij.Because mi and mj originate from the same node, mi and mj mainly depend on the thermal energy flowing through pipes i and j.Therefore, the value of ρij is very small, and the correlation coefficient between mi and mj can be approximated as 0. As can be seen from Lemma 4, mi and mj can be considered to be independent of each other.The flow balance equation for node k can be determined from Figure 2.
Illustration describing the calculation of the mass flow rates correlation coefficient.

 
Combining Equation (11) and Lemma 3 yields the following: Because ρij ≈ 0, Equation ( 12) can be given as: As shown in Appendix B, setting the sum of all thermal loads flowing through pipe i to i   and its variance as  and setting the sum of all thermal loads flowing through pipe j to j   and its variance as yield the following expressions: Because the variance  m of the mass flow through pipe 1 is known from Equation ( 10), the variances of the mass flow rate through pipes adjacent to pipe 1 can be obtained according to Equations ( 14) and ( 15), and the variances of mass flow rates through all other pipes in the network can be obtained in the same way.This process is generalized as follows.
If n pipes are connected to node k and the pipe indices are defined as i1, i2, ..., in, then Equations ( 14) and ( 15) can be established for all mass flow rates in the pipe network as follows: Accordingly, the following equations for pipe i can be obtained: Here, Tstarti and Tendi indicates the temperature at the start and the end of pipe i respectively.Equation ( 19) can be revised according to Equation (20) as follows: which can be rewritten as: If the temperature of node i is Ti, the mass flow rate is from H to node i, and the pipes transmitting the mass flow rates are re-indexed as x1, x2, ..., xk, while the temperatures of the nodes are re-indexed as T .This yields the following for pipe 1 in Figure 1 (i.e., pipe x1): while the following is obtained for the pipe 2 in Figure 1 (i.e., pipe x2): Similarly, this can be extended for an arbitrary pipe xk as follows: Adding Equations ( 23)-( 25) yields the following: which can be written as: Lemma 5. Assume that a continuous random variable X has a probability density function fx(x).It is also assume that a function y = g(x) is monotonous and its inverse function is x = g -1 (x).Accordingly, Y = g(X) is a continuous random variable whose probability density function is: The PDF of a random variable x = m obtained from a normal distribution is: where ui and σi 2 are the mean and variance, respectively.Therefore, the probability density function of y = 1/x = 1/mi is given from Lemma 5 as follows: which can be written in the following form: Based on the form of Equation ( 28), if the mean of y in Equation ( 30) is 1/ui, the standard deviation is (σi•y)/ui, where y = 1/ui.As such, the standard deviation of Equation ( 30) is σi/ui 2 , and Equation (30) can be rewritten as follows: Therefore, if the mean and standard deviation of a random variable x = mi obtained from a normal distribution are, respectively, ui and σi, then y = 1/x = 1/mi approximates a normal distribution, and its mean and standard deviation are 1/ui and σi/ui 2 , respectively.
Similarly, as shown in Appendix C, the correlation coefficient between mass flow rate mk and mi in Figure 2 (ρki) and the correlation coefficient between mass flow rate mk and mj (ρkj) can be obtained as follows: .  32) and (33) and Lemma 6. Accordingly, assuming that the correlation coefficient between mass flow rate passing through pipes x1 and x2 is ρ12, and ρ21 = ρ12, the covariance matrix ∑ of the district heating network can be given as follows: .
When a thermal load fluctuates, the temperature change of the corresponding node is relatively small.Then, the temperatures H T , 27) are desirable for their mean value, where the error is small at this time and can be approximately ignored.

Coupling Elements between Electrical and District Heating Network
The coupling elements acting between electrical and district heating network include cogeneration CHP units, heat pumps, electric boilers, and circulating pumps.Both electrical and thermal energy are supplied by CHP units simultaneously.Heat pumps and electric boilers convert electrical energy into heat.Circulating pumps consume electrical energy to circulate water in the thermal system.These coupling components help to increase the operational flexibility of interconnected electrical-thermal IESs.
CHP units can be divided into two types, depending on whether they employ a fixed thermoelectric ratio (such as gas turbines and reciprocating internal combustion engines) or a variable thermoelectric ratio (such as exhaust steam turbines).A fixed thermoelectric ratio Cm and a variable thermoelectric ratio Cz can be obtained from the electrical energy generation PCHP and the heat generation φCHP of a CHP unit as follows: Here, ηe is the condensation efficiency, and Fin is the fuel input rate of the CHP unit.

Case Study
The 23-node radial district heating network shown in Figure 3 was employed for conducting a case study of the proposed nonlinear analytical algorithm for predicting the probabilistic mass flows of radial district heating network.The CHP source temperature TH is constant at 80°C.The return water temperature To of the load node is constant at 45°C.The ambient temperature Ta is simplified to be a constant, which is set as 10 °C.The remaining parameters of the district heating network are presented in Appendix D.
The mean mass flow rate (μm), standard deviation of the mass flow rate (σm), mean node temperature (μT), and standard deviation of the node temperature (σT) obtained for the test system using the proposed method with those obtained using the Monte Carlo method (simulated 50,000 times) expressed as μm,mcs, σm,mcs, μT,mcs, and σT,mcs, respectively are compared.Assuming that the Monte Carlo simulation values are approximately accurate, the error in our analytical algorithm according to the absolute value differences between the two values obtained

Test 1-Typical Mean and Variances of Network States
1) For part 1 of test 1, the calculated mean and standard deviations of the mass flow rates and node temperatures for selected pipes are shown in Tables 1 and 2, respectively, along with the differences between the values.2) For part 2 of test 1, the calculated mean and standard deviations of the mass flow rates and temperature for selected pipes are shown in Tables 3 and 4, respectively, along with the differences between the values.It can be noted from Tables 1 and 3 that increasing the value of L from 300 m to 1000 m, while holding the mean thermal loads and fluctuations constant, increased both the mean mass flow rates error and the mass flow rates standard deviation error.Nonetheless, the maximum error in the mean mass flow rate was less than 0.03%, while the maximum error in the standard deviation of the mass flow rate was less than 0.004 kg/s.Tables 2 and 4 indicate that similar results were obtained for the mean temperature error and the temperature standard deviation error, where both increased with increasing L, although the maximum mean temperature error was less than 0.002%, while the maximum temperature standard deviation error was less than 0.002 °C.

Test 2-Typical Mean and Variances of Network States
For test 2, the calculated mean and standard deviations of the mass flow rate and temperature for selected pipes are shown in Tables 5 and 6, respectively, along with the differences between the values.In addition, Figures 4(a-d) present the cumulative density functions (CDFs) of mass flow rate through pipe 1 and node 19 temperature under thermal load fluctuations within 10%, 20%, 30%, and 40%, respectively.
It can be noted from Tables 5 and 6 that increasing the thermal load fluctuation with constant L and μ increased the standard deviations of the mass flow rates and the temperature.Nevertheless, for thermal load fluctuations of 10%, 20%, 30%, and 40%, the maximum errors in the mass flow rates standard deviation were 0.0058 kg/s, 0.0127 kg/s, 0.0190 kg/s, and 0.0266 kg/s, respectively, while the maximum errors in the node temperature standard deviation were 0.0028 °C, 0.0072 °C, 0.0138 °C, and 0.0252 °C, respectively.In addition, it can be noted that the mean values of the mass flow rates and the node temperature obtained by the Monte Carlo method decreased with increasing thermal load fluctuation.

Test 3-Probability Density Function of Pipeline Temperature Drops
For test 3, the calculated mean and standard deviations of the pipe temperature drops obtained for selected pipes are shown in Table 7 along with the differences between the values.It can be seen from Table 7 that increasing the thermal load fluctuation with constant L and μ increased the mean pipe temperature drop obtained by Monte Carlo simulations slightly.Meanwhile, increasing the thermal load fluctuation by a factor of 5 with a constant L increased the standard deviation of the pipe temperature drop obtained by Monte Carlo simulations by a factor of 5 generally, and increasing L by a factor of 3 and 1/3 with a constant  increased both the mean and standard deviations of pipe temperature drops by a factor a little less than 3 and 1/3.

Test 4-Control Variable Method
1) From Table 2, It can be noted that the value of σT,mcs is greatest for node 19.Therefore, the value σT,mcs for node 19 with constants μ and σ while varying L is determined, and the results are shown in Figure 5.It can be seen from the figure that the value of σT,mcs for node 19 increases approximately linearly with increasing L. 2) The value σT,mcs for node 19 with constants L and σ while varying μ, is also determined and the results are shown in Figure 6.It can be seen from the figure that the value of σT,mcs for node 19 decreases nonlinearly as μ increases, and approaches zero asymptotically.
. 3) The value σT,mcs for node 19 with constants L and μ while varying σ, is also determined and the results are shown in Figure 7.It can be seen from the figure that the value of σT,mcs for node 19 increases approximately linearly with increasing σ.Nonetheless, the value of σT,mcs for node 19 remains small even at large σ.4) The value σm,mcs for pipe 1 with constant σ while varying both L and μ, is also determined and the results are shown in Figure 8.It can be seen from the figure that the value of L has little influence on the value of σm,mcs for pipe 1, and σm,mcs is essentially unchanged while varying L at any constant value of μ.In contrast, the value of μ has a significant effect on σm,mcs for pipe 1, and σm,mcs increases approximately linearly with increasing μ.5) The value σT,mcs for node 19 with constant σ while varying both L and μ, is also determined and the results are shown in Figure 9.As can be seen from the figure, the value of σT,mcs for node 19 is dependent on both L and μ and increases with increasing L and decreasing μ.Here, L has a small effect on σT,mcs at high μ, but the effect of L is quite significant at low μ, and σT,mcs increases markedly with increasing L.

Model Error Analysis
The values of σm, σm,mcs, and δσ,m for pipe 1 and the values of σT, σT,mcs, and δσ,T for node 19 were compared at a constant value of σ = ±10% while varying both L and μ.
1) When given values of μ with L varied from 100 m to 2000 m in increments of 100 m, 20 sets standard deviation by Monte Carlo and the proposed method respectively can be get, whose average are σm and σm,mcs and the results are shown in Table 8.It can be seen from the table that L has little influence on the standard deviation of pipe 1 mass flow rate under constants μ and σ and that the values of σm and σm,mcs for pipe 1 increase approximately linearly which are shown in Figure 10.  Figure 10.The standard deviations of mass flow rate through pipe 1 using the Monte Carlo method and the proposed analytical method with respect to average thermal load (σ = ±10%).
2) Figure 11 presents the values δσ,T for node 19 with constant σ while varying both L and μ.It can be seen from the figure that the error of the proposed model increases with increasing L and decreasing μ, but the value of δσ,T for node 19 is typically less than 0.0003°C, so the error of the proposed model is generally within an acceptable range.

Figure 11.
Relationship between the differences in the standard deviations obtained for the node temperature of node 19 using the Monte Carlo method and the proposed analytical method with respect to pipe length and average thermal load (σ = ±10%).

Conclusions
This paper developed a nonlinear analytical algorithm for predicting the probabilistic mass flow of a radial district heating network based on the principle of heat transfer and basic pipe network theory.The validity and rationality of the proposed algorithm was verified by application to a 23-node radial district heating network with various pipe lengths under thermal load fluctuations of various magnitudes.The characteristics of the algorithm and the conclusions obtained are given as follows: 1) The proposed algorithm utilizes a nonlinear mass flow model with several reasonable approximations.Consequently, the obtained operating conditions are sufficiently accurate.
2) The algorithm provides probabilistic operational information for district heating network with stochastic heat loads.The algorithm not only provides the variances of the mass flow rate through a pipe network and the node temperatures but also obtains the variances of the pipe temperature drops.3) The computation is efficient because the probabilistic district heating network mass flow model is relatively simple, and our approach does not require repeated nonlinear mass flow calculations.
4) The case study results indicate that the pipe length has little effect on the standard deviations of mass flow rates, while the mean thermal load significantly influences the standard deviations of mass flow rates.
The algorithm proposed in this paper is expected to be very useful for the calculation of district heating network probability and for conducting risk analysis.According to Equations (A5) and (A6), the node temperature changes less when the thermal load fluctuates.Therefore, the following equations can be obtained based on the nature of a normal distribution: .
Accordingly, the following equation can be obtained:  The following equation can be obtained by combining Equations (B3) and (B7):

Appendix C
The following equation can be obtained from the properties of the correlation coefficient: Assuming that ki  is the correlation coefficient between mass flow rate mk and mi, the following equation can be obtained: Because mk = mi + mj and E(X + Y) = E(X) + E(Y), Equation (A15) can be rewritten as follows: (( ) ) ( ) ( ) Because mi and mj are approximately independent, the following equation can be obtained: E(mimj) − E(mi)E(mj) ≈ 0. (A17) Furthermore, from the nature of a normal distribution, the following equation can be obtained: E(mi 2 ) − E(mi) 2 = D(mi).(A18) Finally, combining Equations (A16)-(A18) yields the following: The following equation can be similarly proven: The line parameters of the 23-node district heating network are listed in Table A1, where the thermal load nodes are nodes 7, 8, 10, 11, 12, 14, 15, 16, 19, 20, 21, and

Figure 1 .Lemma 1 . 2 ~
Figure 1.Probabilistic mass flow model of a radial district heating network consisting of a heat source node H, three pipes, and three nodes.Here, the circled values represent pipes and the arrows represent the direction of mass flow rates.

Lemma 7 .
If X = (X1, X2, …, Xn) follows the n-dimensional normal distribution N(a, B) and C is an arbitrary m × n matrix, then Y = C•X follows the m-dimensional normal distribution N(C•a, C•B•C T ), where a and B are the mathematical expectation and covariance matrix of the random variable X, respectively.According to Lemma 7, the probability distribution of k x T (or Ti) in Equation (27) can be obtained from Equations (32)-(34).

Figure 3 .
Figure 3. System diagram of the 23-node heat distribution network employed as a case study.

Figure 8 .
Figure 8. Relationship between the standard deviation of mass flow rate for pipe 1 with respect to pipe length and mean thermal load (Test 4: σ = ±10%).

Figure 9 .
Figure 9. Relationship between the node temperature standard deviation of node 19 with respect to pipe length and mean thermal load (Test 4: σ = ±10%).
g th o f e a ch p ip e (k m ) 0 T h e r m a l lo a d ( M W ) SD of mass flow rate through pipe 1 (kg/s) Difference in SDs for node 19 temperature ( le n g th (k m ) T h e r m a l lo a d ( M W )

Author Contributions:
Methodology, case study & writing original manuscript: G.S. and W.W.; editing & validation: Z.W., H.Z. and S.C.; supervision & review: W.H., Y.W. and Z.Y.Setting the sum of all thermal loads flowing through pipe i to i   in Figure A1, setting the variance of the thermal loads to 2 i    , setting the sum of all thermal loads flowing through pipe j to j   , and setting the variance of the thermal loads to

Table 2 .
Typical mean and standard deviations of node temperature.

Table A1 .
22. Line parameters of the 23-node district heating (H is the thermal source).