Abstract
When applying the widely-used Newton-Raphson (N-R) method to the integrated solution of the hybrid power-water flow problem of the power-water integrated energy system (PW-IES), it is found that there are numerical differences among Jacobian elements, which may yield an ill-conditioned Jacobian matrix and even cause convergence problem. Therefore, this letter proposes to per-unitize the water distribution network (WDN) to reduce the numerical differences and thus to improve the Jacobian matrix conditioning and the integrated N-R solution to the hybrid power-water flow problem.
RECENT years, there are significant efforts worldwide to develop the integrated energy system (IES) or multi-energy system (MES) toward a more efficient, flexible, and decarbonized energy system [
Load flow analysis, also termed as energy flow or multi-energy flow in the field of the IES, is the basic tool for IES planning and operation. A lot of studies have been devoted to the load flow modeling and solution of the IES, including the power-gas IES [
The techniques for solving the load flow of an IES can be generally classified into the decomposed method and the integrated method. In our previous work [
The integrated solution solves all equations of various energy vectors simultaneously and the most widely-used method is the Newton-Raphson (N-R) method [
It is understood that a well-chosen per-unit scheme for power system analysis can simplify the modeling and reduce the computational cost. However, other energy vectors including WDN are usually modeled with actual values. Although per-unitizing the gas and heating networks are stated in [
The main contributions of this letter are suggested as follows: ① the per-unit scheme for the WDN with a focus on base value selection is proposed; ② the per-unitized hydraulic equations and the unified Jacobian matrix for the integrated N-R solution to the hybrid power-water flow problem are developed; and ③ the effects of per-unitized WDN on Jacobian matrix conditioning and the integrated N-R solution are discussed with both theoretical analysis and numerical studies.
The AC power flow equations are used to model the power system as:
(1) |
(2) |
where represents the set of buses connecting to bus i; Vi is the voltage magnitude of bus i; is the angle difference between bus i and bus j; and are the real and imaginary parts of the element in the admittance matrix on bus i and bus j, respectively; and are the active and reactive power generation on bus i, respectively; and and are the active and reactive power loads on bus i, respectively.
A typical WDN consists of reservoirs, tanks, pumps, valves, pipes, and loads. The network can be represented by a set of nodes and edges. The set of nodes consists of the reservoir and tank nodes with known hydraulic heads or flows and other nodes with unknown hydraulic heads. The set of edges represents pipes, pumps, and valves with unknown flows. Valves can be treated as a special kind of pipes. Therefore, the only models of pipes and pumps will be presented below.
The steady-state characteristic of a pipe is mainly described by the hydraulic head loss across the pipe. qpk is the flow rate of pipe k, which connects nodes i and j with hydraulic heads Hi and Hj, respectively. The hydraulic head loss can be expressed as [
(3) |
where the subscript p denotes pipes; n is the variable depending on the hydraulic head loss model; and is the resistance coefficient of pipe k and can be expressed as:
(4) |
where Lpk, Dpk, and Cpk are the length, diameter, and roughness coefficient of pipe k, respectively; and and m are variables depending on the hydraulic head loss model. In this letter, we use the Hazen-Williams equation to describe the hydraulic head loss of pipes, where , m, and n are 10.67, 4.871, and 1.852, respectively [
The widely-used variable-speed centrifugal pump is considered. quk is the flow rate of pump k with the from-node i and the to-node j. The hydraulic head gain of the pump is then modeled as [
(5) |
where the subscript u denotes pumps; , , and are the static head, internal resistance coefficient, and hydraulic head exponent of pump k, respectively; and is the relative speed of pump k, i.e., the ratio of working speed to the nominal speed.
Equations (
(6) |
where qk is the flow rate of pipe k or pump k connecting between node i and node j; and is defined as:
(7) |
Besides (6), the WDN is also governed by the continuity equations, which represent mass conservation requiring that the total flow into a node equals the flow out:
(8) |
where and are the sets of edges (including pipes and pumps) entering and leaving node i, respectively; and is the net flow rate injected into node i.
For a WDN with ne edges connecting to nn unknown-head-nodes and n0 known-head-nodes, we define the incidence matrix Aw relating edges to nodes as:
(9) |
Partitioning into the following two submatrices, we can obtain:
(10) |
where A12 is the matrix relating the edges to the nodes with unknown hydraulic heads; and A10 is the matrix relating the edges to the nodes with known hydraulic heads.
The hydraulic equations of the WDN can then be written in the matrix form as:
(11) |
where q is the vector of edge flow rates; H is the vector of hydraulic heads to be determined; H0 is the vector of known hydraulic heads; qw is the vector of known nodal injected flow rates; and A11 is the diagonal matrix with the elements as:
(12) |
where is the general function of pipe k or pump .
Electricity-driven pumps are the main coupling components between the power system and the WDN. For pump k with and quk, the power consumption Puk is [
(13) |
where is the water density; g is the gravitational acceleration constant; and is the efficiency of pump k.
It should be noted that the power factors of pumps are assumed to be constant and the reactive power consumption Quk of pump k can be determined accordingly. For the pump connecting to node i of the power system, the active and reactive power consumption are incorporated into PLi and in (1) and (2), respectively.
The nominal voltages are usually selected as base voltages, and therefore, the per-unit voltages are near 1.0 under normal operating conditions. However, for the WDN, a wide numerical range of parameters or variables can be found because of the actual values used. For example, the hydraulic heads of the 13-node WDN sample system are found to be in the range of ]m. Other examples are pipe parameters. The pipe lengths and diameters for a typical WDN generally lie in ]m and ]mm [
Similar to the per-unit scheme of the power system, the per-unit value of any quantity of a WDN is defined as the ratio of its actual value to the base value. Therefore, the most important thing for per-unitizing a WDN is the method to select the base values.
To develop the per-unit scheme for WDN analysis, we first make a comparison between the variables of the power system and the WDN, as shown in
Scenario | Parameter | Variable |
---|---|---|
Power system | Voltage | |
Voltage drop | ||
Current | ||
Impedance | ||
Apparent power | ||
WDN | Hydraulic head | |
Head loss | ||
Flow rate | ||
Resistance coefficient | ||
Equivalent power |
To per-unit the power system equations, base power SB and base voltage VB are usually specified, and base current and base impedance are then determined to follow the circuit laws or relations between system variables. To define the per-unit scheme for the WDN analysis, the base flow rate qB and base hydraulic head HB are selected independently, and the base values of resistance coefficient and are determined as:
(14) |
(15) |
where the subscript B denotes the base value of each variable.
The resistance coefficient depends on the individual pump, then the base value for pump k is defined as:
(16) |
With the base values, (11) and (13) can be per-unitized, respectively, as:
(17) |
(18) |
where the superscript denotes the per-unit value; and is defined as:
(19) |
Note that when incorporating pump powers into (1), the base power of Puk should be converted to the base power of the power system , which leads to:
(20) |
Another point worth noting is that the hydraulic equations of a WDN should be per-unitized on a common for all pipes or pumps and a common for all nodes throughout the network.
Similar to the selection of SB and VB of the power system, both qB and can be selected quite arbitrarily. However, a well-chosen and HB can help simplify the calculation and make the per-unit values more informative. In this letter, we suggest to select qB and HB with a reference to the characteristic parameters of pipes and pumps, respectively.
According to [
Diameter (mm) | Flow rate ( |
---|---|
100-400 | 0.005-0.113 |
>400 | 0.113-0.176 |
In the light of
As for HB, the static heads of all the pumps shown in (5) are used by setting ), where m is the total number of pumps.
Equations (
(21) |
where contains the unknown variables of the hybrid power-water flow problem, and θ and V are the vectors of the magnitudes and angles of unknown voltage in the power system, respectively, and q and H are the vectors of flow rates and hydraulic heads of unknowns in the WDN, respectively.
Applying the N-R method to (21), the iteration scheme is:
(22) |
where the superscripts k and denote the current iteration and the iteration; and J is the Jacobian matrix calculated by:
(23) |
where and are the Jacobian matrices of the power system and the WDN, respectively; is the derivative matrix of the power equations with respect to flow rates and hydraulic head; P and Q are the active and reactive power functions, respectively; and D is a diagonal matrix calculated by:
(24) |
Note that Jew in (23) is non-zero due to the power consumption of pumps. More specifically, for bus i supplying the power to the pump k connecting nodes m and n, both PLi in (1) and in (2) are the functions of the unknown hydraulic heads Hm and Hn, and the unknown flow rate qk of the pump.
The N-R iteration in (22) proceeds until the maximum magnitude of is less than the specified tolerance or the iteration count exceeds the maximum number of iterations.
As mentioned above, when expressing hydraulic equations in actual values, there are significant numerical differences between the quantities of the two networks. Such significant numerical differences can also be found in Jacobian elements. More importantly, these differences may yield a great condition number of J, and thus lead to the numerical stability problem and even the convergence problem. Here, we provide a simple illustration of this statement.
is the eigenvalues of J, and the values of satisfy:
(25) |
where E, Ee, and Ew are the identity matrices; and and are the eigenvalues of Jee and Jww, respectively.
Since the proposed per-unit scheme is only used to deal with the hydraulic equations, nothing changes in and therefore only will be discussed below. Jww can be written in (26) and are the eigenvalues of , B, and A, respectively, where .
(26) |
According to the min-max theorem and interlacing theorem of eigenvalues [
(27) |
Considering that ① D is a positive and diagonal matrix, which means and ; and ② , and , which means τn is a small negative number while τ1 is a small positive number. We can conclude from (27) that a significant numerical difference in the elements of D will yield a great range of . The 2-norm condition number of Jww, denoted by , is proven to be . One can decrease by reducing the numerical differences among , which is the motivation of this letter to per-unitize the WDN.
Expressing the hydraulic equations in per-unit as in (17) and (18), the hybrid power-water flow problem in per-unit and its N-R iteration scheme, which have similar forms as (21) and (22), can be obtained. The Jacobian matrix in per-unit can be written as:
(28) |
where
(29) |
In the case studies, we will verify numerically that the condition number of the Jacobian matrix can be reduced significantly with the proposed per-unit scheme.
For the N-R method, the initialization should be done carefully. A common practice for initializing power flow is the flat start, i.e., set all unknown voltage and the unknown angle . This flat start can be extended to a per-unitized WDN by setting the per-unit value of all unknown flow rates and . Due to the special consideration in selecting qB and HB as presented above, this initialization is usually proven to be a good start.
Two test systems are studied below to illustrate the advantage of the proposed per-unit scheme in the analysis of hybrid power-water flow. The first system is composed of the modified IEEE 33-bus system incorporating two photovoltaic sources and a 13-node WDN, which is referred to as IEEE 33-W13 below. The second system, denoted by IEEE 33-W27, is composed of the same IEEE 33-bus system and a 27-node WDN.
In
Value | IEEE 33-W13 | IEEE 33-W27 | ||
---|---|---|---|---|
Max | Min | Max | Min | |
Actual value | 14632.304 | -21.494 | 576.046 | -21.494 |
Per-unit value | 21.494 | -21.494 | 21.494 | -21.494 |

Fig. 1 Colormaps of Jacobian matrix of IEEE 33-W27 system. (a) Actual value. (b) Per-unit value.
In
Both
Value | IEEE 33-W13 | IEEE 33-W27 |
---|---|---|
Actual value |
4.4932×1 |
4.7478×1 |
Per-unit value |
2.4641×1 |
2.4641×1 |
For hydraulic analysis, a large variation in pipe lengths, diameters, and roughness may result in an ill-conditioned system [
Case 1: the original parameters of the WDN are considered.
Case 2: the length of pipe W1-W2 increases to three times the original length.
Case 3: the diameter of pipe W1-W2 reduces to one third of the original diameter.
Case | Actual value | Per-unit value | ||||
---|---|---|---|---|---|---|
к(J) | k | T (ms) | к( | k | T (ms) | |
1 |
4.7478×1 | 5 | 7 |
2.4641×1 | 5 | 6.5 |
2 |
5.3035×1 | 5 | 7 |
2.4641×1 | 5 | 7.0 |
3 |
9.8541×1 | - | - |
6.9513×1 | 6 | 8.1 |
Note: “-” represents non-convergent.
In
Case 3 in actual value is the only unsolved case. The eigenvalue distributions of the Jacobian matrix of Case 3 are shown in

Fig. 2 Eigenvalue distribution. (a) J. (b)
Both
This letter proposes a per-unit scheme for the WDN to improve the integrated N-R solution for the hybrid power-water flow problem of the PW-IES. The effects of per-unitizing on the Jacobian matrix conditioning are analyzed theoretically and the improvements in N-R solution are verified numerically through two test systems.
References
P. Mancarella, “MES (multi-energy systems): an overview of concepts and evaluation models,” Energy, vol. 65, no. 2, pp. 1-17, Jan. 2014. [Baidu Scholar]
Q. Zeng, J. Fang, J. Li et al., “Steady-state analysis of the integrated natural gas and electric power system with bi-directional energy conversion,” Applied Energy, vol. 184, pp. 1483-1492, Dec. 2016. [Baidu Scholar]
Q. Sun, Q. Dong, S. You et al., “A unified energy flow analysis considering initial guesses in complex multi-energy carrier systems,” Energy, vol. 213, p. 118812, Dec. 2020. [Baidu Scholar]
Y. Chen, J. Zhao and J. Ma, “Fast decoupled multi-energy flow calculation for integrated energy system, ” Journal of Modern Power Systems and Clean Energy, vol. 8, no. 5, pp. 951-960, Sept. 2020. [Baidu Scholar]
X. Li, X. Zhao, and L. Yang, “Combined power-water flow analysis of regional integrated electricity and water networks,” in Proceedings of 2018 2nd IEEE Conference on Energy Internet and Energy System Integration, Beijing, China, Oct. 2018, pp. 1-6. [Baidu Scholar]
E. Abraham and I. Stoianov, “Constraint-preconditioned inexact newton method for hydraulic simulation of large-scale water distribution networks,” IEEE Transactions on Control of Network Systems, vol. 4, no. 3, pp. 610-619, Sept. 2017. [Baidu Scholar]
X. Yan, “Water conveyance and distribution engineering,” in Water Supply Engineering, 4nd Ed. Beijing: China Architecture & Building Press, 1999, pp. 36-38. [Baidu Scholar]
F. Zhang. (2011, May). Matrix theory. [Online]. Available: https://link.springer.com/book/10.1007/978-1-4614-1099-7. [Baidu Scholar]