Pipe Hydraulic Resistances Identification of District Heating Networks Based on Matrix Analysis

District heating networks (DHNs) are essential municipal infrastructure in the north of China. Obtaining accurate resistance characteristics is a critical step to improve the operating regulation level of DHNs. In this paper, pipe hydraulic resistances (PHRs) are introduced to express the resistance characteristics. A hydraulic model of a DHN can be established by using observed data of pressures and discharges. The boundary nodes are taken as observed sites. After establishing a matrix equation and analyzing the rank of its coefficient matrix, the authors propose a method to determine all the PHRs uniquely, by using a small number of observed sites and operating conditions. Furthermore, when observed errors are introduced, the adverse impact can be weakened by increasing the number of operating conditions and the accuracy of observed devices. When the observed error ranges are 1% and 0.5%, the results show that the average relative errors of identified PHRs are 2.4% and 1.1% respectively, which can be acceptable in engineering. Then, a loop DHN can be transformed into several branch DHNs, which are identified individually.


Introduction
District heating systems are necessary facilities to ensure indoor environment quality in severe cold and cold regions, which cover almost 70% of the area in China. By 2016, the total building floor area with heating facilities had exceeded 7 billion m 2 [1]. District heating networks (DHNs) undertake the task of delivering hot water in the systems. Since the heat transfer performance could be affected by flow conditions [2], it is very useful to obtain accurate resistance characteristics and hydraulic conditions of DHNs. In this paper, the authors focus on hydraulic conditions of DHNs, but not thermal conditions. The establishment of accurate resistance models of DHNs is the foundation of smart heating, energy efficiency, and enhancing the quality of building interior environments. Varying from water distribution networks, the water in DHNs is in a circulation state. Since there is large and continuous water flow, DHNs are sensitive to resistance characteristics. Thus, it is more crucial to identify resistance characteristics in DHNs.
As there is no mathematical expression difference between DHNs and water distribution networks, the identification of resistance characteristics of DHNs can utilize the research results of water distribution networks [3]. Identifying resistance characteristics of water distribution networks usually needs observed data of pressures and flows when the topology is determined. Liu et al. express the resistance characteristics in a water distribution network by pipe hydraulic resistances (PHRs), being related to pipe flows and node pressures directly [4]. The PHRs have been defined by Energies 2020, 13 Wang et al. [5], which are related to Darcy friction factor, length, internal diameter, the loss coefficients of local losses, and the density of water. The accuracy of PHRs identification results will be affected by the condition of observed sites [6,7]. However, few researches focused on the arrangement of observed sites from observability analysis. Limited by the number of observed sites, the current researches need to solve underdetermined systems during PHRs identification. A unique solution to the PHRs identification problem cannot be determined directly. To decrease the number of independent variables and obtain overdetermined systems, many methods have been introduced. Scholars employed a global gradient algorithm to solve governing equations, in which the number of unknowns could be reduced by grouping the pipes according to similar Hazen-Williams coefficients [8,9]. Besides not reflecting the local resistance, the validity of the grouping assumption cannot be guaranteed. To obtain more observed data, multiple operating conditions (OCs) are often needed [10]. Even so, identifying PHRs may still be underdetermined, and the unique solution cannot be obtained theoretically. To express PHRs uniquely, constraints are often added to narrow the ranges of variables [11,12]. However, it is difficult to set the ranges accurately. Liu et al. [13] identified PHRs without the abovementioned simplification. This method built overdetermined systems by using observed data under multiple OCs. Since a large number of observed sites must be needed, it is difficult to apply this method in actual DHNs. In addition, observed errors are also important factors in PHRs identification. Optimization methods based on the least square principle are often used to reduce the influence of observed errors [14,15]. However, the optimization model still needs to solve underdetermined systems. Optimization methods are excessively time-consuming, and are not applicable to large-scale networks, generally.
There are a few researches on PHRs identification of DHNs. Zhou et al. [16] determined PHRs uniquely by using all the node pressures and some pipe flows. Since a large number of observed sites are needed, this method is not applicable to the observed conditions of the current DHNs. Wang et al. [5] identified PHRs of a branch DHN (without loops) by using node pressures and discharges under a large number of OCs. Since the quantitative relationship between the number of independent equations and the number of OCs is not involved, the achievement of unique PHRs will not be guaranteed, and the number of OCs required is relatively large.
Identifying PHRs for both branch and loop DHNs is worth studying, based on the actual observed conditions and operating features. To reduce the number of OCs required and improve the accuracy of the PHRs identification results, a PHRs identification method based on matrix analysis is proposed, laying the foundation for building more accurate resistance models of district heating systems.

Methods
A DHN is a closed network, composed of supply pipelines, return pipelines, and boundary nodes. The boundary nodes include heat sources and heating substations. The PHRs identification of a DHN usually treats the supply pipelines and the return pipelines independently. Since there is no theoretical difference between the supply and return pipelines, the authors take the supply pipelines as an example to analyze. Considering the current observed conditions, pressures and discharges of boundary nodes will be employed to identify PHRs of a DHN.

Basic Equations
In this paper, the governing equations of a branch DHN can be expressed by the mass continuity and energy conservation equations. Actually, the mass continuity equation is also the mass conservation equation, expressing the quantitative relationships between pipe flows and node discharges (being also called node outlet flows) [5].
where A is the basic incidence matrix; A = (a ij ) n×n , (i = 1,2, . . . , n), an element a ij can be described as follows: the values −1 and +1 indicate that the flow in the pipeline is directed toward and away from Energies 2020, 13, 3007 3 of 11 the node, respectively; otherwise, the value is 0 [4]. G is the vector of pipe flows; and Q is the vector of node discharges. A DHN consists of internal nodes and boundary nodes. An internal node is a connecting node, and the corresponding node discharge equals zero. A boundary node represents a heat source or a heating substation, and the corresponding node discharge equals the flow in the heat source or the heating substation. Generally, the discharges of boundary nodes are available. In a branch DHN, when a heat source node is selected as the reference node, the nonsingularity of the basic incidence matrix can be guaranteed [4]. Therefore, the pipe flows can be solved by Equation (1) when the discharges of all the boundary nodes are available. Next, the energy conservation equation is employed to express the quantitative relationship among the node pressures, pipe flows, and PHRs [5].
where P is the vector of node pressures; G diag is the diagonal matrix of pipe flows; and S is the vector of PHRs. To find the pressures of internal nodes involved in Equation (2), the vector of node pressures and the basic incidence matrix need to be expressed as partitioned form, The subscripts b and i represent the variables of boundary nodes and internal nodes, respectively. Then, Equation (2) can be partitioned according to different node types, and expressed as In DHNs, the velocity of the internal flow is almost 0.5-3 m/s. The flow condition of the water is turbulent flow. Thus, the PHRs are independent of the velocity and flow, and the PHRs will not change with variable OCs [17]. To eliminate dependent variables in the matrix equations, the authors establish the quantitative relationships of pipe flows between the kth and the first OCs.
where k is the order number of an OC; C is a diagonal matrix, and a diagonal element equals the flow's square ratio of the corresponding pipe. Suppose that there are N OCs, N partitioned energy conservation equations shown in Equation (3) can be expressed. Then, by using Equation (4), pipe flows and PHRs can be replaced by (N − 1) diagonal matrices C. Combined with the above N equations in partitioned matrix form, a matrix equation with pressures of internal nodes as independent variables can be written as Equation (5).
Since the pressures of the boundary nodes can be observed, the coefficient matrix and non-homogeneous term in Equation (5)  pseudoinverse to Equation (6). (Actually, the values of PHRs can also be solved by using the pipe flows and node pressures under any OC.) . . .

Solution Analysis
Solving pressures of internal nodes is a key part to identify PHRs. Thus, the authors discuss the solution of Equation (5).

Moore-Penrose Pseudoinverse of Matrix
For D · x = b, a matrix D + is the Moore-Penrose pseudoinverse of D (a non-zero matrix), which has the same dimensions as D T . If D is full-column rank, D · x = b will have a unique solution.
The Moore-Penrose pseudoinverse of D can be expressed as The solution of Moore-Penrose pseudoinverse can be expressed as x = D + · b, which satisfies the norm D · x − b equals zero [18]. (If D is a nonsingular matrix, the matrix D + is also the inverse matrix of D.)

Rank Analysis of Matrix
Whether the coefficient matrix in Equation (5) is full-column rank will depend on the dependence of its column vectors. If the coefficient matrix is regarded as a partitioned matrix, a submatrix composed of columns two to N can be considered as a partitioned diagonal matrix. Apparently, every column vector in the submatrix cannot be expressed by other column vectors linearly. A diagonal matrix −C(k) (k = 1, 2, . . . , N − 1) left multiplying A T i is equivalent to row transformation, multiplying different scalars (the principal diagonal elements of the matrix −C(k)) to the corresponding rows of A T i . Since A T i is composed of columns corresponding to every internal node, there must be more than one element equaling 1 or −1 in every column of A T i . Thus, there must be more than one element not equaling zero in any column of −C(k)·A T i . When the independent OCs can be guaranteed, the scalars to different rows are not equal to each other. Thus, the column one in the partitioned coefficient matrix must not be linearly expressed by the submatrix (composed of columns two to N in the coefficient matrix). In the partitioned coefficient matrix, appending the column one to the submatrix will increase its rank. In this case, the coefficient matrix in Equation (5) is full-column rank. Thus, Equation (5) has a unique solution. The pressures of internal nodes under every OC can be determined uniquely.

Quantity Requirement of Independent OCs
If more OCs can be provided, the number of algebraic equations in the matrix equation (Equation (5)) will increase. The coefficient matrix of Equation (5) may be full-column rank. If the coefficient matrix of Equation (5) is full-column rank, unique solutions can be calculated directly. Even when the observed errors cannot be neglected, Moore-Penrose pseudoinverse solutions of Equation (5) can be introduced to compute the unique solutions [18]. A large number of actual engineering cases show that these unique solutions can generally satisfy engineering requirements [4] (in the following part, the case study results also verify the effectiveness of this method).
Suppose there are n b boundary nodes (heat sources and heating substations) and n i internal nodes in a branch DHN. If the number of nodes or pipes is n, then n b + n i = n and n b − n i ≥ 1 (in a branch DHN, besides a reference node, the number of nodes equals the number of pipes). To ensure the coefficient matrix of Equation (5) is full-column rank, the minimum quantity of independent OCs N needs to satisfy (N − 1) · n ≥ N · n i . Then, the minimum integer of N is two. Thus, for a branch DHN, pressures of internal nodes can be determined uniquely by pressures and discharges of boundary Energies 2020, 13, 3007 5 of 11 nodes under two independent OCs. Then, all the PHRs can be solved by the calculation results of pipe flows and node pressures.

PHRs Identification of a Loop DHN
In a loop DHN, the number of pipes is greater than the number of nodes (except the reference node). In a basic incidence matrix corresponding to a loop DHN, the number of columns is greater than the number of rows. Thus, the basic incidence matrix must be singular, and the pipe flows cannot be reached by known node discharges. If observed devices are equipped additionally in some internal nodes, the proposed method may still be applied. However, this requirement is beyond the current technical conditions of DHNs. As the values of PHRs are independent of the velocity and flow (caused by different layouts or OCs of DHNs), a transforming method is introduced. To realize PHRs identification under the current technical conditions, a method for transforming a loop DHN into several branch DHNs is introduced. By shutting off different pipes (valves) of one loop in the original loop DHN, several different branch DHNs appear. PHRs of these branch DHNs can be identified by the proposed method individually. When the PHRs identified cover all the pipes, the PHRs identification of the loop DHN will be achieved.

Case Study
In this section, two numerical examples will be used to verify the feasibility of the proposed method. The authors first take a branch DHN as an example to introduce the identification processes. Furthermore, the authors also realize the PHRs identification of a loop DHN.

A Branch DHN
An example branch DHN contains 12 nodes (including one reference node, n0) and 11 pipes. The topology of the DHN is shown in Figure 1.
According to the previous definition, a basic incidence matrix (A) corresponding to the DHN (supply pipelines) shown in Figure 1 can be written as follows.
Pipes information in the DHN is listed in Table 1. By using the setting values of PHRs and flows, the pressures and discharges of the boundary nodes, being also called operating data and assumed as known "observed values" in the following context, can be computed according to hydraulic calculation equations. Those pressures and discharges are assumed as available "observed values".
PHRs identification processes are expressed by a numerical example without considering observed errors. In this situation, the PHRs will be identified when the known operating data are assumed to be accurate. Then, the influence of observed data with errors are studied in accordance with the observed conditions of actual DHNs.
Energies 2020, 13, 3007 6 of 11 method. The authors first take a branch DHN as an example to introduce the identification processes. Furthermore, the authors also realize the PHRs identification of a loop DHN.

A Branch DHN
An example branch DHN contains 12 nodes (including one reference node, n0) and 11 pipes. The topology of the DHN is shown in Figure 1. According to the previous definition, a basic incidence matrix (A) corresponding to the DHN (supply pipelines) shown in Figure 1 can be written as follows.  Table 1. By using the setting values of PHRs and flows, the pressures and discharges of the boundary nodes, being also called operating data and assumed as known "observed values" in the following context, can be computed according to hydraulic calculation equations. Those pressures and discharges are assumed as available "observed values".

PHRs Identification Processes
For the branch DHN shown in Figure 1, given the operating data under two independent OCs, the theoretical unique values of PHRs can be identified. The operating data (generated by hydraulic calculation) are employed to replace actual DHN data, and are listed in Table 2 (the relative pressure at the reference node n0 is 110 m). The PHRs identification is realized in Matlab environment, using calculation tools (such as functions ABS, DIAG, INV, and PINV). The major procedures are listed below.
(1) Calculate the pipe flows under the two OCs by using the discharges shown in Table 2 (Equation (1)).
(3) Calculate the rank of the coefficient matrix in Equation (5). When the coefficient matrix is full-column rank, the theoretical unique solutions of Equation (5) can be solved. The pressures of internal nodes can be found in the solution of Moore-Penrose pseudoinverse.
The identification results can verify that the coefficient matrix in Equation (5) is full-column rank, and can be guaranteed when operating data under two OCs are available. The identification results of PHRs are completely equal to the setting values of PHRs shown in Table 1, illustrating the validity of this PHRs identification method.

Analysis of Observed Errors
This section analyzes PHRs identification when the observed errors cannot be neglected. It is assumed that all the observed errors obey the normal distribution. The two error ranges (1% and 0.5%) are selected, respectively. A Matlab function, Normrnd, is called to generate random numbers which obey the normal distribution. Then, the random numbers generated are added to numerical calculation values (without errors) of operating data (pressures and discharges). The sum of the two parts expresses observed data with errors, which can be used as known operating data during the PHRs identification processes. Those operating data with errors are shown in Tables 3 and 4. In the processes of PHRs identification, it makes no difference whether there are observed errors. Theoretically, any two or more independent OCs (such as OC 1-2) can guarantee to determine all the PHRs values uniquely. However, due to the observed errors, another OC (OC 3) is added. In this section, the operating data under OC 1-2 and OC 1-3 are analyzed, respectively.
To evaluate the accuracy of identified PHRs, a relative error of an identified PHR can be expressed as where δ is the relative error of an identified PHR; s is the identified value of a PHR; and s 0 is the setting value of a PHR. The relative errors of PHRs identification related to different error ranges (1% and 0.5%) are shown in Figures 2 and 3, respectively.  When the observed error range is 1%, the relative errors of PHRs are shown in Figure 2. The maximum and average relative errors under two OCs are 41.4% and 11.2% respectively, and the maximum and average relative errors under three OCs are 5.5% and 2.4% respectively. Similarly, when the observed error range is 0.5%, the relative errors of PHRs are shown in Figure 3. The maximum and average relative errors under two OCs are 24.4% and 8.1% respectively, and the maximum and average relative errors under three OCs are 4.9% and 1.1% respectively. When only two OCs are selected to determine all the PHRs, if the setting values of PHRs are relatively smaller (also with less resistance), the relative errors of these PHRs are usually higher (such as p5 and p6 in Figure 2, and p6, p7, and p8 in Figure 3). Figures 2 and 3 indicate that the relative errors will decrease as the number of OCs increases. By comparing between Figures 2 and 3, it can be concluded that the relative errors will decrease as the accuracy of observed devices increases (the error ranges decrease).  When the observed error range is 1%, the relative errors of PHRs are shown in Figure 2. The maximum and average relative errors under two OCs are 41.4% and 11.2% respectively, and the maximum and average relative errors under three OCs are 5.5% and 2.4% respectively. Similarly, when the observed error range is 0.5%, the relative errors of PHRs are shown in Figure 3. The maximum and average relative errors under two OCs are 24.4% and 8.1% respectively, and the maximum and average relative errors under three OCs are 4.9% and 1.1% respectively. When only two OCs are selected to determine all the PHRs, if the setting values of PHRs are relatively smaller (also with less resistance), the relative errors of these PHRs are usually higher (such as p5 and p6 in Figure 2, and p6, p7, and p8 in Figure 3). Figures 2 and 3 indicate that the relative errors will decrease as the number of OCs increases. By comparing between Figures 2 and 3, it can be concluded that the relative errors will decrease as the accuracy of observed devices increases (the error ranges decrease). When the observed error range is 1%, the relative errors of PHRs are shown in Figure 2. The maximum and average relative errors under two OCs are 41.4% and 11.2% respectively, and the maximum and average relative errors under three OCs are 5.5% and 2.4% respectively. Similarly, when the observed error range is 0.5%, the relative errors of PHRs are shown in Figure 3. The maximum and average relative errors under two OCs are 24.4% and 8.1% respectively, and the maximum and average relative errors under three OCs are 4.9% and 1.1% respectively. When only two OCs are selected to determine all the PHRs, if the setting values of PHRs are relatively smaller (also with less resistance), the relative errors of these PHRs are usually higher (such as p5 and p6 in Figure 2, and p6, p7, and p8 in Figure 3). Figures 2 and 3 indicate that the relative errors will decrease as the number of OCs increases. By comparing between Figures 2 and 3, it can be concluded that the relative errors will decrease as the accuracy of observed devices increases (the error ranges decrease).

A Loop DHN
An example loop DHN contains 14 nodes (including one reference node, n0) and 14 pipes, as listed in Table 5. The topology of the example loop DHN is shown in Figure 4.

A Loop DHN
An example loop DHN contains 14 nodes (including one reference node, n0) and 14 pipes, as listed in Table 5. The topology of the example loop DHN is shown in Figure 4.  In this example loop DHN, p2 and p10 are shut off, respectively. By shutting off p2, a branch network (BN 1) is determined and shown in Figure 5a. In BN 1, PHRs of p5, p6, p7, p8, p9, p10, p11, p12, and p13 are identified. Similarly, by shutting off p10, another branch network (BN 2) is determined, and shown in Figure 5b. In BN 2, PHRs of p1, p2, p3, p4, p5, p8, p9, p13, and p14 are identified. Actually, for the identification processes, there is no difference between different shutting schemes (for instance, if p4 and p12 are shut off, the other different two branch networks will appear. PHRs of the two branch DHNs will be identified similarly). In this example loop DHN, p2 and p10 are shut off, respectively. By shutting off p2, a branch network (BN 1) is determined and shown in Figure 5a. In BN 1, PHRs of p5, p6, p7, p8, p9, p10, p11, p12, and p13 are identified. Similarly, by shutting off p10, another branch network (BN 2) is determined, and shown in Figure 5b. In BN 2, PHRs of p1, p2, p3, p4, p5, p8, p9, p13, and p14 are identified. Actually, for the identification processes, there is no difference between different shutting schemes (for instance, if p4 and p12 are shut off, the other different two branch networks will appear. PHRs of the two branch DHNs will be identified similarly).
Identifying PHRs of the example loop DHN shown in Figure 4 is equivalent to identifying PHRs of the two branch DHNs shown in Figure 5. Identifying every PHR can be achieved by using operating data under two independent OCs for both BN 1 and BN 2. The identification processes are in accordance with the procedures applied in Section 3.1. The identification results of PHRs are listed in columns BN 1 and BN 2 in Table 5. Furthermore, the authors choose mean values to express the identification results of PHRs shown in column BN 1-2 in Table 5 (in the identification processes, p15, p16, p17, and p18 are virtual pipelines, which are not included in the original DHN. Thus, the authors will not list those identification results). The maximum and average relative errors of PHRs identification results compared with the setting values, shown in the last column of Table 5, are 5.8% and 0.9%, respectively. The result shows that the proposed method for a loop DHN is feasible, and has high accuracy. This method can realize PHRs identification by operating data under multiple OCs, and there is no need to equip extra observed devices on the internal nodes in loop DHNs. Identifying PHRs of the example loop DHN shown in Figure 4 is equivalent to identifying PHRs of the two branch DHNs shown in Figure 5. Identifying every PHR can be achieved by using operating data under two independent OCs for both BN 1 and BN 2. The identification processes are in accordance with the procedures applied in Section 3.1. The identification results of PHRs are listed in columns BN 1 and BN 2 in Table 5. Furthermore, the authors choose mean values to express the identification results of PHRs shown in column BN 1-2 in Table 5 (in the identification processes, p15, p16, p17, and p18 are virtual pipelines, which are not included in the original DHN. Thus, the authors will not list those identification results).
The maximum and average relative errors of PHRs identification results compared with the setting values, shown in the last column of Table 5, are 5.8% and 0.9%, respectively. The result shows that the proposed method for a loop DHN is feasible, and has high accuracy. This method can realize PHRs identification by operating data under multiple OCs, and there is no need to equip extra observed devices on the internal nodes in loop DHNs.

Conclusions
This paper first establishes a PHRs identification method for a branch DHN. This method can be realized only by using pressures and discharges of boundary nodes. By rank analysis of the matrix equation, the minimum number of OCs making the PHRs identification problem an overdetermined system is determined. When the coefficient matrix of the energy conservation equation corresponding to all the OCs is full-column rank, there is a unique solution to the matrix equation. The pressures of internal nodes under every OC can be found in the solution of Moore-Penrose pseudoinverse. Then, the values of PHRs can be obtained. Theoretically, the results of PHRs are the real solution of the problem. For a loop DHN, the PHRs identification can be realized by identifying several branch DHNs. The above research results have been realized by numerical examples. The following conclusions have been drawn: (1) For a branch DHN, all the PHRs can be determined uniquely by pressures and discharges of boundary nodes under two independent OCs.
(2) As the number of OCs and the observed device accuracies increase, the relative errors of PHRs identification results will decrease. When the error range is 1%, the average relative errors of PHRs identification results under two and three OCs are 11.2% and 2.4%, respectively. When the error range is 0.5%, the average relative errors of PHRs identification results under two and three OCs are 8.1% and 1.1%, respectively.
(3) For a loop DHN, all the PHRs can still be determined uniquely by shutting off different pipes and building several different branch DHNs. The average relative error of PHRs identification

Conclusions
This paper first establishes a PHRs identification method for a branch DHN. This method can be realized only by using pressures and discharges of boundary nodes. By rank analysis of the matrix equation, the minimum number of OCs making the PHRs identification problem an overdetermined system is determined. When the coefficient matrix of the energy conservation equation corresponding to all the OCs is full-column rank, there is a unique solution to the matrix equation. The pressures of internal nodes under every OC can be found in the solution of Moore-Penrose pseudoinverse. Then, the values of PHRs can be obtained. Theoretically, the results of PHRs are the real solution of the problem. For a loop DHN, the PHRs identification can be realized by identifying several branch DHNs. The above research results have been realized by numerical examples. The following conclusions have been drawn: (1) For a branch DHN, all the PHRs can be determined uniquely by pressures and discharges of boundary nodes under two independent OCs. (2) As the number of OCs and the observed device accuracies increase, the relative errors of PHRs identification results will decrease. When the error range is 1%, the average relative errors of PHRs identification results under two and three OCs are 11.2% and 2.4%, respectively. When the error range is 0.5%, the average relative errors of PHRs identification results under two and three OCs are 8.1% and 1.1%, respectively. (3) For a loop DHN, all the PHRs can still be determined uniquely by shutting off different pipes and building several different branch DHNs. The average relative error of PHRs identification results is 0.9%, which can illustrate the feasibility and validity of the method.
In addition, this PHRs identification method of DHNs is based on the governing equations (mass continuity and energy conservation equations). Thus, the fundamental principle of the proposed method can also be applied to other fluid supply networks, such as gas and water distribution networks.
Author Contributions: Conceptualization: Y.L. and P.L.; funding acquisition: Y.L. and P.L.; methodology: Y.L. and P.W.; writing-original draft: Y.L.; writing-review and editing: Y.L. and P.W. All authors have read and approved the final version of the manuscript.