State Estimation of District Heating Networks PDF - Applied Thermal Engineering

Applied Thermal Engineering 73 (2014) 1211e1221 Contents lists available at ScienceDirect Applied Thermal Engineering journal homepage: www.e...

Applied Thermal Engineering 73 (2014) 1211e1221 Contents lists available at ScienceDirect Applied Thermal Engineering journal homepage: www.elsevier.com/locate/apthermeng State estimation of district heating network based on customer measurements Tingting Fang*, Risto Lahdelma Department of Energy Technology, School of Engineering, Aalto University, Otakaari 4, Espoo 02150, Finland h i g h l i g h t s
Estimating the state of a district heating network based on customer measurements.
The new model can be adapted to treat extra measurements to improve robustness.
State estimation can help with various operational network managements.
Volume flows of both tree-structured and cyclic networks are discussed.
Automated remote meter reading facilitates better accuracy for the model. a r t i c l e i n f o a b s t r a c t Article history: District heating (DH) has been widely used in many European countries since the beginning of 20th Received 27 November 2013 century. In Finland approximately half of the heating market was covered by DH by 2011. For better use Accepted 1 September 2014 of district heating, real-time monitoring of the water flows, temperatures and heat losses can be an Available online 16 September 2014 effective vehicle to manage DH networks. Automated remote meter reading is a new technology that was applied since 2005 and will cover the Helsinki region within 10 years. Automated meter reading enables Keywords: computing the heat losses more accurately than before. In this paper, a new model is designed to es- District heating timate the water flows, temperatures, and heat losses in different parts of the network using automated State estimation Automated meter reading hourly meter readings for customers. Assuming a connected graph topology for the network, the flow Least squares estimation equations form a determined and potentially nonlinear system that we solve iteratively. Then, based on Heat losses the computed water flows and customer temperature measurements, we express the water temperature and heat loss in each pipe as an over-determined linear system and solve it in the least squares (LSQ) sense which ensures the solution with the smallest 2-norm error vector. We illustrate the method with a real life DH network based on hourly temperature and flow measurements for one week. The model can be applied to arbitrary DH networks. © 2014 Elsevier Ltd. All rights reserved. 1. Introduction Approximately 74% of DH is produced in combined heat and power (CHP) plants. DH technology features high energy efficiency of In 1877, the first district heating (DH) network was established the CHP production, low level of emissions, customer-friendliness, in Lockport (New York, USA). Since the beginning of the 20th possibility to use various fuels, and utilization of heat storages. DH century, DH has been widely utilized in many European Countries can have a renewable energy supply stock, and also allows for such as Germany, Denmark, Holland, Belgium, Sweden and Finland. regulating emissions in centralized plants. The high energy effi- The first DH network of Finland was built in 1940 in Helsinki. In ciency of DH systems in comparison to other heating options has Finland, DH has covered almost half of the heating market by 2011 been shown in various studies. For example and describe the. In most big cities, the market share has hit 90% [4,5]. general advantages of DH (and cooling) systems. Ref. discusses the role of DH in reducing CO2 and NOx emissions. Refs. [9e11] discuss DH under different special conditions, such as together with renewable energy, in sparsely populated areas, and with * Corresponding author. geothermal source. Surplus heat from industrial processes can also E-mail addresses: tingting.fang@aalto.fi, [email protected] (T. Fang), risto.lahdelma@aalto.fi (R. Lahdelma). be utilized in DH systems, then DH enables further enhancement of http://dx.doi.org/10.1016/j.applthermaleng.2014.09.003 1359-4311/© 2014 Elsevier Ltd. All rights reserved. 1212 T. Fang, R. Lahdelma / Applied Thermal Engineering 73 (2014) 1211e1221 Nomenclature Tis supply water temperature at node i ( C) (measured for i 2 Nmea) Tir return water temperature at node i ( C) (measured for Abbreviations i 2 Nmea) DH district heating s TðiÞj temperature of incoming supply water from node i to LSQ least squares (estimation) branching node j ( C) (j 2 Nbranch) r TðjÞi temperature of incoming return water from node j to Index sets branching node i ( C) (i 2 Nbranch) N set of nodes in the DH network N ¼ Nmea ∪ Ncomp Nmea set of measured nodes Parameters Ncomp set of computed nodes r density of water (kg/m3) Nbranch set of nodes where three or more branches meet C specific heat capacity of water (J/(kg  C)) A set of DH pipes connecting nodes di measured demand (supply) of water at the node i O set of cycles (m3/s) li,j the length of pipe (m) Variables Si,j combined heat loss factor for the pipes between nodes v_ i;j volume flow rate in supply and return pipes i and j (W/m  C) (m3/s) T g ground temperature ( C) Fsi;j heat loss for the supply pipe (W/m) ( 2 A) fi,j Darcy friction factor for pipe (1) Fri;j heat loss for the return pipe (W/m) efficient resource usage. For the customer, DH is an economical and management of DH. The same is true also for district cooling as well user-friendly heating solution. However, DH also has some chal- as for gas and electricity distribution. Many companies are lenges, such as large investment costs, as well as heat losses of the moving towards remote meter reading, which means that the network during transportation of the heat. collected data of the consumed heat are sent to the energy com- In Finland, the DH network is a dual pipe system where heat is pany automatically. In the past, it was possible to measure only the delivered to customers as heated water. The supply pipe takes the annual heat losses of the overall network since the automated hot water to the customers and the return pipe brings the cooled meter reading measurements were not available. With remote water back to the heating plant. The water is circulated with pumps meter reading, it is possible to collect from the customers a large which are located at the heating plant; for some larger networks amount of more accurate energy consumption data in real time. they also stay along the DH network. Usually the network's con- This source of information enables the calculation of heat losses struction pressure is 1.6 MPa. The customers are connected to the more accurately than before. Since 2005, remote meter reading network with indirect connections, which means that the heat has been applied in some parts of Helsinki region. This technology from the DH system is transferred through heat exchangers in is expected to be applied throughout Helsinki within 10 years. In buildings to secondary radiator circuits and also to heat up do- the short term, remote reading increases the cost of DH network, mestic water. but in the long run costs will be saved. At the moment, remote Heat loss in a DH system is mainly caused by conduction of heat meter readings are used mainly for billing purposes but not from the pipes into the ground. A small fraction of the heat is also monitored consecutively. In this paper, we take into consideration conducted from the supply pipe into the return pipe. Inter-pipe the automated remote meter reading data from customers, which heat conduction does not cause heat loss, but reduces the heat enable us to compute the hourly heat losses for each pipe. Remote transfer capacity of the network a little. Many studies have been meter reading can also be applied by energy companies for done in Refs. [12e15] to evaluate the heat losses and analyse the detecting faults in the network and identifying inefficient heat efficiency for DH systems. Kitanovski and Poredos calculated usage of customers. An important use of measurement data is to heat losses in the heat distribution network as 8e10%. The heat loss forecast the demand more accurately. can be as large as 10e20% in small networks with an average pipe Existing DH network calculation software systems do not utilize (inner) diameter of 50 mm. In large networks with average pipe the abundance of measurement data to compute an optimal esti- diameter of 150 mm the heat loss is around 4e10%. Larger pipe mate of the DH network state. In this paper we design a new model surface area in comparison to transferred heat will cause the larger to compute the optimal least squares estimate for the state of the losses in small networks. It's impractical to manually measure the DH network based on hourly customer consumption measure- heat losses of each pipe in a large network. Therefore, a new ments. This means that the state estimation model will output the computational method to estimate the heat losses in different parts volume flow rate and heat loss in each pipe along with the supply of the network using all available measurements is desirable. and return temperature at each node based on the input mea- Three main measurements, volume flow rate of the DH water as surements of volume flow, and supply & return water temperatures well as the temperatures of the supply and return water, are used of the customers. With this output it is possible to monitor the by the heat consumption meter to compute the heat consumption. network state more accurately and in real-time. The state infor- Earlier, the customers sent their heat meter readings to the energy mation is useful for managing and operating the network more company on a regular basis (e.g. once a year). One disadvantage of efficiently. More detailed temporal and spatial network state in- this approach was that the readings were always incomplete. In formation can also help in planning the DH production optimally. addition, the distribution of consumed heat energy could not be An earlier version of state estimation for a small cycle-free network resolved. Installing automated remote meter reading for all cus- was presented in Ref.. This study extends the model for net- tomers is a recent trend in DH technology. In the near future, EU works with cycles and demonstrates the method using a part of a legislation will require more advanced measurement and real-life network of a Finnish municipal DH company.  T. Fang, R. Lahdelma / Applied Thermal Engineering 73 (2014) 1211e1221 1213 This paper is organized as follows. A real life DH network is should equal the demand, water supply d1 can be computed by introduced in Section 2. Section 3 illustrates how to determine the summing all the customer demands. From Fig. 2 we can see that the volume flow and temperature as well as heat loss. Computational demand variation of different customers is to some extent coinci- results for the real life network are presented in Section 4. Section 5 dent due to space heating depending on outdoor temperature and contains the conclusions and ideas for future work. partly non-coincident due to random use of the hot tap water. From Fig. 3 we can see that due to heat consumption by the customers, the return water is about 40  C cooler than the supply water. 2. Network representation In addition, the known quantities also consist of the pipe data including length, diameter, insulation thickness, thermal conduc- Some definitions regarding district heating networks should be tivity of the insulation and thermal properties of the ground. The clarified before we introduce the network: pipe data was provided by the DH company. Definition 1. A node is the junction of two or more branches. 3. Estimation based on customer measurements Definition 2. A branch is an element of the network connecting two nodes. One branch corresponds to two pipes (supply and re- The state estimation is implemented by the following three turn pipes). steps: Definition 3. A branching node is the junction of three or more branches. 1. Identification of the topology of the network to verify that the network is connected and also to locate any cycles. This means Definition 4. A intermediate node is the junction of any two or identifying the set of nodes N (partitioned into measured nodes more than two branches. Branching nodes are subset of interme- Nmea & nodes to be computed Ncomp) and the set of pipes A as diate nodes. well as the set of cycles O. Definition 5. A cycle is a closed circuit of consecutive branches. 2. Estimation of the volume flow rate v_ i;j for each pipe based on the volume flow measurements at measured nodes as well as Definition 6. A customer node is typically a leaf node in the related pipe data (e.g. pipe lengths, inner diameter, heat trans- network. mission factors) and flow properties (e.g. Reynolds numbers, We assume that the DH network forms a connected graph. friction factors etc.). In practice a DH network is mostly tree-like, but it may contain a 3. Estimation of the temperatures at each computed node (Tis small number of cycles. Cycles make the network more robust, andTir ), and the heat loss (Fsi;j and Fri;j ) for each pipe based on because in case of pipe failure an alternative route to distribute heat measured temperatures as well as the volume flows computed may be available. To illustrate the model, Fig. 1 shows a real life DH in step 2. network for a suburb in Helsinki region. Node 1 represents the heat plant (or entry point to the network). Intermediate nodes 20, 22, 23, The computations in steps 2 and 3 have been implemented in 25, 26 and 29 are located along a cycle which is marked with green Matlab. dashed arrows. The network contains 14 customer nodes 4, 12, 15 etc. that are located as leaf nodes in the network. Fig. 1 also shows 3.1. Estimation of water flows the nominal water flow directions in the supply pipes. In reality, the water direction can sometimes be reversed, in particular in pipes We first consider the estimation of water flows in a network along a cycle. Also, if the network contains multiple supply points, without cycles and then with cycles. the flow directions may depend on the supply and demand of water. 3.1.1. Volume flow in a tree-like network Input data consists of hourly measurements for the water de- For a tree-structured network without loss of water, the volume mands and temperatures at customer nodes for one week (168 h), flow in each supply pipe is identical to the corresponding return which are shown in Figs. 2 and 3, respectively. As the water supply pipe (but opposite in direction). This means computing the flows Fig. 1. A real life DH network with 14 customers. 1214 T. Fang, R. Lahdelma / Applied Thermal Engineering 73 (2014) 1211e1221 Fig. 2. Water demand measurements for customer nodes of the sample real life network. for the supply network is sufficient to determine all the volume The same solution can be obtained without matrix computa- flows. tions by summing up recursively the water demand starting from For each node in the network, an equation can be written to the customer nodes along the branches towards the supply node. state that the incoming water subtracted by the outgoing water If also the flow measurement for the supply node or any inter- must equal the (measured) water consumption or supply at that mediate nodes are available, (2) is an over-determined system. node: An over-determined system has more equations than variables. The X X extra equation can be either redundant or inconsistent. If the v_ j;i  v_ i;j ¼ di c i2N (1) supply measurement differs from the overall consumption, the < j;i > 2R < i;j > 2R system is inconsistent and over-determined. Such difference could be caused by water leak or (more typically) by measurement errors. At customer nodes, there will be a positive (non-negative) de- The inconsistent over-determined system has no exact solution, but mand di. At intermediate (non-customer) nodes, the consumption it can be solved in the least squares (LSQ) sense. is zero. The water supply at the power plant can be considered as a negative consumption. We can rewrite (1) in matrix form as !  T 1 T ! v_ ¼ F F F d: (4) ! ! F v_ ¼ d ; (2) The LSQ solution is the optimal estimate that minimizes the Euclidian norm of the error vector. Matlab will solve system (2) where F is the connectivity matrix. The connectivity matrix is a automatically in the LSQ sense when the system is over- matrix where each row represents a node of the network and each determined. Existing DH calculation software do not handle such column represents a branch (with a nominal direction) of the over-determined systems. network. The matrix will contain þ1 where the node is the desti- nation of the branch, 1 where the node is the origin of the branch 3.1.2. Volume flow in a network with cycles and 0 elsewhere. A DH network may contain a small number of cycles. Here we When the graph is a tree (connected and non-cyclic) with n focus on traditional DH networks with cycles. Special ring topology nodes, system (2) contains n equations and n-1 variables. One of the networks are out of the scope of this paper. Mathematical model- equations is redundant, because the supply at node 1 is the sum of ling ring topology with mass flow control is discussed in Refs. all the demands. Omitting the equation corresponding to the sup- [21,22]. Cycles make the flow system more complicated to solve, as ply node (or any one of the equations) yields a fully determined water can take multiple paths between nodes, and the volume linear equation system that has a unique solution flows in the supply and return pipes are not necessarily identical. In ! ! this case solving the volume flows requires flow dynamic v_ ¼ F 1 d : (3) modelling.  T. Fang, R. Lahdelma / Applied Thermal Engineering 73 (2014) 1211e1221 1215 Fig. 3. Supply and return temperature measurements for customer nodes of the sample real life network. A cyclic network can be formed from a tree-like network by ! X Pð v_ Þ≡ ±Dpi;j ðv_ i;j Þ ¼ 0; (5) adding branches between the nodes of the tree-like network. In a < i;j > 2Ok tree-like network, every pair of nodes is connected by a unique path. Each branch added to the tree creates an additional path where Dpi,j is the pressure drop or gain for the branch along the between some of the nodes. A connected cyclic network can also be cycle. The sign of the pressure drop Dpi,j is chosen based on the flow reduced back into a tree by removing cycle-forming branches. direction in the corresponding pipe. The pressure drops in the flow Each cyclic branch added to a tree-like network adds one new direction and increases against the flow. The pressure drop in each variable to the system (2), but does not add new equations. As a pipe Dpi;j ðv_ i;j Þ is a non-linear function of the volume flow rate (see result, the number of excess variables in system (2) equals the Appendix). By combining the pressure balance equations for cycles number of cycles. This means that system (2) becomes under- with the linear flow Equation (2) gives the general volume flow determined, and cannot alone be used to solve all flows in the system model as network. An under-determined system has fewer independent 8 ! ! equations than variables. The system will be under- < F v_ ¼ d  ! X determined both in presence and absence of the nth measure- : Pk ð v_ Þ ¼ ±Dpi;j v_ i;j ¼ 0; k ¼ 1; …; M ; (6) ment of the supply node. Consider as an example the cycle in the < i;j > 2Ok sample network between nodes 20, 22, 25, 29, 26, 23. We can select an arbitrary branch along the cycle and eliminate the correspond- where M is the number of cycles in the network. For our sample ing variable from (2) by setting the flow to a constant value. The network, M ¼ 1, and the linear under-determined system contains remaining flows can then be solved from the system (2) as in the 50 equations and 51 variables. case of a tree-like network. In the general case with multiple cycles, The flows can be solved iteratively from (6) as described in the the flow variable of each cycle-forming branch can be set to an next section. The pressure in a pipe depends quadratically on the arbitrary constant value and the remaining flows can be solved volume flow rate, and linearly on the Darcy friction factor fi,j (see from (2). Appendix). However, the friction factor depends on the volume flow Additional constraints are needed to uniquely solve all flows in a in a complex non-linear manner. With large values for Reynolds cyclic network. A cycle is a closed circuit of consecutive braches. For number (turbulent flow) the friction factor is nearly constant, which each cycle Ok, we can write a pressure balance equation, which means that the pressure drop is close to quadratic function of the flow. states that the sum of the pressure drops and gains in the pipes In a network with cycles, the volume flow in the supply and along the cycle must be zero: return networks are not necessarily identical. This is caused by 1216 T. Fang, R. Lahdelma / Applied Thermal Engineering 73 (2014) 1211e1221 different water temperature in the supply and return pipes (the where m and n is now the number of equations and unknowns difference is about 40  C), which affects the viscosity, density, respectively. Formula (7) is actually the one dimension case of Reynolds number and friction factor. Normally, the flow difference Formula (8). The initial value of v_ can be chosen similar to the one- between the supply and return systems is small, and it is sufficient dimensional case. to estimate the flows only for the supply system and assume the same (but reverse) flow in the return system. For more accurate results, it is possible to solve the flows in the supply and return systems separately using different temperature data solved by the 3.2. Estimation of temperatures and heat losses temperature & heat loss model which is discussed in Section 3.2. Heat loss from DH pipes is proportional to the temperature difference between the surrounding ground and the water in the 3.1.3. Solving the volume flow system for a network with cycles pipes. The water temperature will drop exponentially over time. As With multiple cycles, we can solve the volume flow system the water travels at constant speed in the pipe, the temperature will using the iterative NewtoneRaphson method. In the single cycle also decrease exponentially as a function of the travelled distance. case, Newton's method can be applied. Newton's method is the one The blue curve in Fig. 4 shows that for a supply pipe. The flows and dimensional special case of the NewtoneRaphson method. their directions have been determined earlier, and from here on we The sample network contains one cycle. To solve the system, we assume that the direction is from i towards j in the supply system apply Newton's iterative method. The iterative formula which can and opposite in the return system. Because the temperature drops be found in Ref. is written as very little along a single pipe segment between two nodes , we   s in the pipe as P v_ ðkÞ can approximate the average water temperature Ti;j v_ ðkþ1Þ ¼ v_ ðkÞ  0  ðkÞ  _ (7) the average of the temperature at the endpoints. The heat loss Fsi;j of P v_ a pipe can be approximated to be proportional to the length of the pipe and also proportional to the difference between the average Here, P ($) is the single pressure balance equation in (5). The iter- temperature Ti;js and the ground temperature Tg. Next we form the ation stops when the difference in the flow variable between two temperature and heat loss equations based on the above consecutive iterations is smaller than a predetermined value, i.e. approximations.  ðkþ1Þ  v_  v_ ðkÞ  < d. The single cycle model contains n nodes and n variables, n-1 linear flow equations and one non-linear pressure balance equa- tion. We choose a flow variable corresponding to an arbitrary 3.2.1. Temperature equations branch along the cycle as the variable to iterate in (6). With a fixed We first consider a small element in a supply pipe. The value for this variable, the remaining n-1 variables can be solved time for travelling the pipe length is from the n-1 linear constraints (as discussed earlier). To start the iteration an arbitrary value for v_ ð0Þ can be selected. For example, li;j v_ ð0Þ ¼ 0 is possible. The iteration will converge faster, if the starting ti;j ¼ : (10) ai;j value is close to the actual value. One strategy to achieve this is to select as the iterated branch one where the flow may be close to where ai,j is the velocity of the flow in the pipe. The energy balance zero. That would be a narrow pipe far from the supply node. for this element of length dl can be written as Another strategy to speed up iteration is to assume an initial value ð0Þ for the friction factors fi;j along the cycle and then to solve the   li;j initial volume flow of the selected pipe from a resulting quadratic dlDrC Tis  Tjs ¼ dlFsi;j ; (11) ai;j equation. The initial friction factors can be computed based on average water temperature and intermediate value for Reynolds where D is the cross sectional area of pipe. Rearranging the equa- number. Then, at current iteration, the friction factors are recom- tion and substituting v_ i;j ¼ Dai,j we obtain a linear dependency puted based on the flows computed at the previous iteration. between the temperature difference at the end nodes of the pipe With multiple cycles, the system (6) contains multiple pressure and the heat loss of the pipe balance equations PðvÞ. _ The argument v_ is a vector of flow variables chosen so that they correspond to one branch of each cycle v. _ Then we use the NewtoneRaphson method to solve the volume flows. The iterative NewtoneRaphson method is based on the equation !ðkþ1Þ !ðkÞ !ðkÞ  !ðkÞ  v_ ¼ v_  J 1 v_ P v_ ; (8)   where J v_ / ðkÞ is the Jacobian matrix which contains the partial !ðkÞ derivatives of the function Pð v_ Þ and denoted as below 2 3 vP1 vP1 6 vv_ /   6 1 vv_ n 7 7 6 7 J v_ / ðkÞ ¼ 6 6 « 1 « 7 7; (9) 6 7 4 vPm vPm 5 / vv_ 1 vv_ n Fig. 4. Water temperature as function of travelled distance (blue) and linear approx- imation (red). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)  T. Fang, R. Lahdelma / Applied Thermal Engineering 73 (2014) 1211e1221 1217 ground and between the supply and return pipes. Ti;j s þ T r =2 Fsi;j li;j i;j Tis  Tjs ¼ c < i; j > 2A (12) stands for the average water temperature in the two pipes. Tg is the v_ i;j rC ground temperature ranging from approximately 5  C in the southern part to 2  C in the northern part of Finland. We break When node j is the destination of multiple flows from different down the above equation for the supply pipe and return system branches, each branch brings water of different temperature into respectively: the node. This means that we cannot use the same variable Tjs in   (12) to represent the end-temperature of each pipe. Instead, we Fsi;j ¼ Si;j Ti;j s  Tg ; (18) must introduce pipe-specific end-node temperature variables TðiÞj s. Fig. 5 illustrates this in the case of two incoming flows from nodes   im and in into node j for both the supply and return system. Fri;j ¼ Si;j Ti;j r  Tg : (19) To consider multiple incoming flows towards node i, Equation (12) becomes A part of the heat loss from the supply pipe is recovered in the return pipe as the supply and return pipes are normally installed Fsi;j li;j Tis  TðiÞj s ¼ c < i; j > 2A (13) close to each other. Therefore a slightly smaller S should be used for v_ i;j rC the return pipe, but for simplicity we disregard this. Then the water temperature in each pipe can be computed as an average of the s ¼ T s , when j is not a target of multiple flows. At nodes where TðiÞj j temperature at the end nodes to obtain again linear dependency target for multiple incoming flows, the temperature of the mixed between node temperatures and heat losses. water is the weighted average of the water temperatures in the As with temperature equations, TðjÞis and T r must be used to ðiÞj incoming pipes: replace Tis and Tjr when multiple flows meet at a node. When a node X v_ i;j is not a target of multiple flows, we have TðjÞi s ¼ T s and T r ¼ T r. Tjs ¼ P s TðiÞj : (14) i ðiÞj j ij < i;j > 2R v_ i;j Therefore, (18) and (19) can be written as ij < i;j > 2R s þ Ts ! Similarly, the temperature equations for return system can be TðjÞi j expressed as Fsi;j ¼ Si;j T g c < i; j > 2A; (20) 2 Fri;j li;j Tjr  TðrjÞi ¼ ; c < i; j > A (15) ! v_ j;i rC Tir þ TðiÞj r Fri;j ¼ Si;j T g c < i; j > 2A (21) r ¼ T r , when node j is not a target of multiple flows. Also, the TðjÞi 2 i temperature of mixed water in the return system is computed similarly: The combined heat conduction coefficient Si,j can be calculated X based on the resistance for heat conduction in the pipe insulation, v_ j;i Tir ¼ P r TðjÞi : (16) the resistance between the pipes, and in the ground as described in jj < i;j > 2R jj < i;j > 2R v_ j;i Ref.. In the long run (years to decades) the insulation may lose its properties, but in the short term it heat conductivity will remain constant. The heat conductivity of the ground is subject to random 3.2.2. Heat loss equations variations, mainly due to changes in ground moisture. To calculate heat losses in DH pipes, analytical methods and explicit solutions for the regular cases can be found in Ref.. Heat loss from a pair of insulated buried DH pipes is generally calculated as s þ Tr ! Ti;j i;j g Fi;j ¼ 2Si;j T ; (17) 2 where the factor 2 accounts for the two pipes. Si,j is the combined heat loss factor for the pipes between the ith and the jth node, which takes into account heat conductivity of the insulation, the Fig. 6. Structure of the temperature-heat loss estimation model, H1, H2, H3, H4, H5, H6, Fig. 5. Temperature variables around a branching node (a trivial case for both supply right hand side vector and the size of this system are determined by the and return system). measurements. 1218 T. Fang, R. Lahdelma / Applied Thermal Engineering 73 (2014) 1211e1221 8 > > Fsi;j li;j > > Tis  TðsiÞj ¼ c < i; j > 2A > > v_ i;j rC > > > > > > > > Fri;j li;j > > Tjr  TðrjÞi ¼ ; c < i; j > A > > v_ j;i rC > > > > ! > > > > TðsjÞi þ Tjs > > Fsi;j ¼ Si;j  Tg c < i; j > 2A > > 2 > > > < ! Tir þ TðriÞj (22) > > F r ¼ S  T g c < i; j > 2A > > i;j i;j 2 > > > > > > X v_ > > > > T s ¼ Xi;j TðsiÞj cj2ftarget of multiple supply flowsg > > j v_ i;j > > ij < i;j > 2R > > > > ij < i;j > 2R > > > > X v_ > > Tir ¼ Xj;i TðrjÞi ci2ftarget of multiple return flowsg > > > > v_ > : jj < i;j > 2R j;i jj < i;j > 2R 3.2.3. Solving the temperature heat loss system with actual measurements of flows, temperatures and heat losses To summarize, each branch of the DH network contributes 4 at different pipes. However, such extra measurements are not linear equations. Two of these are temperature equations; two are available for real-life networks. Instead, we have simulated a heat losses equations. In addition, target nodes of multiple segment of the real life network using Grades Heating version 2.57 incoming flows contribute to linear equations. These equations are which manages detailed calculation for the DH network in terms of for supply and return pipes respectively: temperatures, volume flows, mass flows, heat losses, pressure We rewrite the above linear system in matrix form as drops, etc. Grades Heating was originally developed by Process ! ! Vision Ltd, which was merged with Enoro in 2012. Grades Hx ¼ z: (23) Heating has been used widely in Nordic countries, for example in After rearranging the model (22), the coefficient matrix of the modelling the Oulu DH system (1600 GWh). Because Grades linear system can be written as a partitioned matrix which is Heating cannot calculate an over-determined system, we entered to shown in Fig. 6. Grades heating as input the supply temperature of the heating The temperature-heat loss system (23) contains 2 temperature plant computed by our model, and the temperature differences and variables per unknown node, one additional temperature variable heat loads measured at customer nodes. After simulation we per branching pipe, and 2 heat loss variables per branch. In the case compared the flows, heat losses and temperatures calculated by of a tree-like network (connected and non-cyclic) with all Grades Heating with our model. The results were found sufficiently customer nodes measured, the supply-side of the system is at least consistent to validate our model. Small differences were explained determined, and normally over-determined if at least one branch- by the fact that our model assumes constant water density based on ing node exists. However, the return-side is only determined. If also average water temperature in the network. To obtain identical re- the supply node or some other intermediate nodes are measured, sults, we would have to iterate solving the flow system with then both the supply and return systems are over-determined. For a computed temperatures and resolving the heat losses and tem- cyclic network which can be formed from a tree-like network by peratures based on the new flows. adding branches between the nodes of the tree-like network, the temperature and heat loss model (23) will be over-determined, 4. Computational results because each additional cycle adds more equations than variables to the model. Our sample DH network has 218 equations and 204 As there are 102 pipes and 51 nodes in the sample real life unknown variables (103 equations and 90 unknowns for supply network, it is not possible to show the results for all the pipes and system; 115 equations and 114 unknowns for return system). nodes in this paper. The computational results for one customer The LSQ solution of the over-determined linear system (23) is path including pipes , , and were  1 selected to show in Figs. 7 and 8. ! ! x ¼ HT H HT z : (24) Fig. 7 shows both supply and return temperatures for nodes 20, 22, 24, 27 and 31 (node 31 is measured). The supply water tem- We form the coefficient matrix and right hand side vector pro- perature remains rather steady, as it is controlled by the heat plant. grammatically based on the network's structure data. The supply In contrast, the return temperature varies more, as it depends on and return systems can be solved independently. Solving two the variable heat consumption of the customers. The return tem- smaller systems can be more efficient than solving one large perature curves for nodes 24, 27 and 31 show a similar pattern as system. the flow in the supply system doesn't branch after node 22. There are different ways to validate the model. To verify that the Monitoring the temperature in different parts of the network is implementation is correct, we have tested the model on several valuable to guarantee proper functioning. We can see that the artificial networks. Ideally, we would like to apply the model on supply temperature at each node is high enough to deliver heat to several real-life networks and compare the computational results each customer substation. The return temperatures are around  T. Fang, R. Lahdelma / Applied Thermal Engineering 73 (2014) 1211e1221 1219 Fig. 7. Estimated supply and return water temperatures for nodes in a customer path from a real life DH network. 40  C which is ideal. Too high return temperature would indicate 5. Conclusion and discussion operating problems such as cooling errors at substations, short circuit flows without control by thermal valves, or too low supply In this paper, an estimation model was developed to estimate temperatures at substations because of high heat loss. Too rapid the water flows, temperatures and heat losses in a DH network temperature changes over time steps will result in thermal based on customer measurements. expansion and contraction of the pipes causing them to wear out The developed model, which manages detailed calculation of faster. Also, the temperature difference at waters mixing at branch water flows, temperatures and heat losses in DH pipes, enables nodes should not be too large, as that causes wear-out in the pipe robust and accurate DH network state estimation. The estimation joints. model allows tracking the state of the DH network much more The estimated heat losses and volume flows for all 4 þ 4 pipes accurately than before. DH companies can apply state estimation to are shown in Fig. 8. The heat losses are reasonably low, which in- rapidly aid in identifying and locating various problems in the dicates that the network is in good shape. The heat losses from the network or in the measurements. Reliable and automatic DH return pipes are significantly smaller than from the supply pipes. network state estimation is a prerequisite for various operational This is consistent with the lower temperature of return water. Also, network management and optimization tasks. The state estimation the return pipes express much more variation than the supply pipes model can assist optimal production planning in a variety of as they are directly affected by customer behaviour. A significant different situations. For example, the heat company may want to variation or trend in the heat losses or volume flow may indicate optimize how much heat they should produce at different heat problems either in the network or in the measurements. In plants to best satisfy the demand in different parts of the DH particular, large inconsistency residue error in the flow system (2) network at the minimum cost. Such optimization could be imple- or (6) could indicate a leakage in the network. mented on top of the state estimation model. 1220 T. Fang, R. Lahdelma / Applied Thermal Engineering 73 (2014) 1211e1221 Fig. 8. Estimated heat losses, volume flows for all the pipes in a customer path from the sample real life DH network. We illustrated the model using a real life network with 51 nodes, 102 pipes and one cycle, but it can be used for arbitrary networks. 8li;j rv_ i;j 2 However, real life networks may result in quite large estimation Dpi;j ðv_ i;j Þ ¼ fi;j : (A.1) dIi;j 5p2 models. With larger networks, it may be necessary to implement special sparsity exploiting techniques to solve the estimation here dIi;j is the inner diameter of the pipe , fi,j is the Darcy models. To improve the accuracy of the model, we could iterate friction factor for pipe , r is the density of water at the corre- between the volume flow and temperature & heat loss models. The sponding temperature, and li,j is the pipe length. estimated temperatures from the temperature & heat loss model Apart from volume flow v_ i;j, the only variable to be determined can be used to re-compute the water density, viscosity, Reynolds in the above equation is the friction factor fi,j. The Moody diagram number and volume flows. The more accurate volume flows can presents how the friction factor depends on the Reynolds number, again be fed into the temperature & heat loss model for better pipe diameter and roughness in both turbulent and laminar flow. In accuracy. practice, flows in DH pipes are mostly turbulent. For turbulent flow, An advantage with the DH state estimation model is that it can the friction factor can be solved from the implicit equation for treat extra measurements at intermediate nodes, making the model turbulent flow below. more over-determined and robust with respect to measurement errors. One direction for future research is to develop methods to 0 1 εi;j assess the accuracy of the measurements and estimated DH 1 B dI 2:51 C qffiffiffiffiffi ¼ 2 log10 @ i;j þ qffiffiffiffiffiA; (A.2) network state based on extra measurements. 3:7 fi;j Rei;j fi;j where εi,j is the roughness of pipe (0.04 mm for the steel Acknowledgements pipes), and Re is the Reynolds number. This equation can be solved uniquely by binary search. The Reynolds number is defined as This work is supported by companies participating in the RICES- project. Special thanks to Helsingin Energia (the municipal energy 4v_ i;j company of Helsinki) for providing the data. Thanks are also given Rei;j ¼ ; (A.3) pdIi;j vi;j to the research milieu at Aalto University. where vi,j is the kinematic viscosity for the flow in pipe. The viscosity and density of DH water depends somewhat on the water Appendix temperature. Assuming an average temperature for the DH and using corresponding values for r and vi,j in Equations (A.1) and (A.3) To determine the pressure drop in each pipe, an analytical is an acceptable approximation. More accurate calculations require function can be found in Ref. , which is written as information about water temperature in different pipes, but this  T. Fang, R. Lahdelma / Applied Thermal Engineering 73 (2014) 1211e1221 1221 information can be obtained only later after we solve the heat loss K. Çomaklı, B. Yüksel, O. € Çomaklı, Evaluation of energy and exergy losses in district heating network, Appl. Therm. Eng. 24 (7) (2004) 1009e1017. and temperature system. A. Poredos, A. Kitanovski, Exergy loss as a basis for the price of thermal en- ergy, Energy Convers. Manag. 43 (2002) 2163e2173. References A. Dalla Rosa, H. Li, S. Svendsen, Method for optimal design of pipes for low- energy district heating, with focus on heat losses, Energy 36 (5) (2011) D.K. Baker, S.A. Sherif, Heat transfer optimization of a district heating system 2407e2418. using search methods, Int. J. Energy Res. 21 (1997) 233e252. M.A. Rosen, M.N. Le, I. Dincer, Efficiency analysis of a cogeneration and district Energiateollisuus ry (Finnish Energy Industries), District Heating Handbook energy system, Appl. Therm. Eng. 25 (1) (2005) 147e159. (in Finnish: Kaukola€mmo € n k€ asikirja), Kirjapaino Libris Oy (printed by Libris), European Union, Energy Efficiency Directive, Access method: ec.europa.eu/ Helsinki, 2006. energy/efficiency/eed/eed_en.htm (referred 14.10.2013). District Heating and Cooling Country by Country e 2005 Survey, Euroheat & K. Kontu, T. Fang, R. Lahdelma, Forecasting district heating consumption based Power, Brussels, 2005. on customer measurements, in: Proc. IASTED Conf. on Power and Energy Official Statistics of Finland (OSF): Energy consumption (e-publication), Systems and Applications, Las Vegas, USA, 2012, pp. 210e216. ISSN¼1798-6869, Helsinki, Statistics Finland (referred: 31.7.2012), Access T. Fang, K. Kontu, R. Lahdelma, Estimation of the state of a district heating method: www.stat.fi/til/ekul/index_en.html. network based on consumption measurements, in: Proc. IASTED Conf. on Helsingin Energia (Helsinki Energy), Access method: www.helen.fi (referred Power and Energy Systems and Applications, Las Vegas, USA, 2012, pp. 6e13. 14.10.2013). G. Chartrand, Introduction to Graph Theory, McGraw-Hill Education Pvt B. Rezaie, M.A. Rosen, District heating and cooling: review of technology and Limited, India, 2006. potential enhancements, Appl. Energy 93 (2012) 2e10. G.H. Golub, C.F. Van Loan, Matrix Computations, fourth ed., JHU Press, S. Frederiksen, S. Werner, District Heating and Cooling, first ed., Exaktap- Maryland, 2013. M. Kuosa, K. Kontu, T. Ma €kila €, M. Lampinen, R. Lahdelma, Static study of rinting AB, Sweden, 2013. M.E. Ossebard, A.J.M. Van Wijk, M.T. Van Wees, Heat supply in the traditional and ring networks and the use of mass flow control in district Netherlands: a systems analysis of costs. Exergy efficiency, CO2 and NOx heating applications, Appl. Therm. Eng. 54 (2) (2013) 450e459. M. Kuosa, M. Aalto, M. El Haj Assad, T. Ma €kila €, M. Lampinen, R. Lahdelma, emissions, Energy 22 (1997) 1087e1098. I.B. Hassine, U. Eicker, Impact of load structure variation and solar thermal Study of a district heating system with the ring network technology and plate energy integration on an existing district heating network, Appl. Therm. Eng. heat exchanger in a consumer substation, Energy Build. 80 (2014) 276e289. 50 (2013) 1437e1446. F.P. Incropera, D.P. Dewitt, Fundamentals of Heat and Mass Transfer, John E. Bowitz, M. Dang Trong, The social cost of district heating in a sparsely- Wiley, New York, 1990. populated country, Energy Policy 29 (2001) 1163e1173. J. Stoer, R. Bulirsch, Introduction to Numerical Analysis, third ed., 2010. I. Yabanova, A. Kecebas, Development of ANN model for geothermal district German. heating system and a novel PID-based control strategy, Appl. Therm. Eng. 51 Grades Heating, Calculation System for District Heating Distribution Net- (2013) 908e916. works, Enoro website (visited 31.08.2014): http://www.enoro.com/.