Journal of Modern Power Systems and Clean Energy

ISSN 2196-5625 CN 32-1884/TK

网刊加载中。。。

使用Chrome浏览器效果最佳,继续浏览,你可能不会看到最佳的展示效果,

确定继续浏览么?

复制成功,请在其他浏览器进行阅读

Per-unit Scheme for Integrated Solution of Hybrid Power-water Flow Problem  PDF

  • Xia Zhao
  • Hong Tan
  • Luo Wang
  • Mingyi Sun
the State Key Laboratory of Power Transmission Equipment & System Security and New Technology, Chongqing University, Chongqing, China

Updated:2023-03-24

DOI:10.35833/MPCE.2021.000453

  • Full Text
  • Figs & Tabs
  • References
  • Authors
  • About
CITE
OUTLINE

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.

I. Introduction

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 [

1]. The IES seeks potential synergies among electricity, natural gas, heating, and cooling networks, and therefore suggests to plan and operate various energy vectors from a more integrated perspective. In this context, there is a growing interest in the interdependences between the power and water as well as the coordination of the power system and water distribution network (WDN).

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 [

2], the power-heating IES, and the power-gas-heating IES [3], [4].

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 [

5], the decomposed solution to the hybrid power-water flow problem of the power-water IES (PW-IES) is proposed and this work focuses on the integrated solution.

The integrated solution solves all equations of various energy vectors simultaneously and the most widely-used method is the Newton-Raphson (N-R) method [

2]-[3], where the unified Jacobian matrix is developed and updated during the iterations. Reference [4] proposes a fast decoupled multi-energy flow calculation method, where the Jacobian matrix is replaced by a diagonal and constant Jacobian matrix. Many studies have demonstrated the good applicability of the N-R method to the load flow solution of the IES. However, when applying the N-R method to the PW-IES, numerical differences among Jacobian elements are found because the power equations are expressed in per-unit value while the hydraulic equations of the WDN are in actual value. These numerical differences may yield an ill-conditioned Jacobian matrix, which may result in the numerical stability problem and even the convergence problem of the N-R solution. Therefore, we propose to per-unitize the WDN to improve the integrated N-R solution of the hybrid power-water problem.

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 [

2] and [3], respectively, the base value selection in per-unitizing and the effects of per-unit modeling on the Jacobian matrix as well as the N-R solutions are not fully addressed.

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.

II. Hybrid Power-water Flow Modeling of PW-IES

A. Power Flow Equation

The AC power flow equations are used to model the power system as:

PGi-PLi=VijiVjGijcos θij+Bijsin θij (1)
QGi-QLi=VijiVjGijsin θij-Bijcos θij (2)

where ji represents the set of buses connecting to bus i; Vi is the voltage magnitude of bus i; θij is the angle difference between bus i and bus j; Gij and Bij are the real and imaginary parts of the element in the admittance matrix on bus i and bus j, respectively; PGi and QGi are the active and reactive power generation on bus i, respectively; and PLi and QLi are the active and reactive power loads on bus i, respectively.

B. Hydraulic Equation of WDN

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 ΔHij can be expressed as [

6]:

ΔHij=Hi-Hj=γpkqpkqpkn-1 (3)

where the subscript p denotes pipes; n is the variable depending on the hydraulic head loss model; and γpk is the resistance coefficient of pipe k and can be expressed as:

γpk=αLpk/CpknDpkm (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 [

6].

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]:

ΔHuk=Hj-Hi=ωuk2H0,uk-rukωuk2-mukqukmuk (5)

where the subscript u denotes pumps; H0,uk, ruk, and muk are the static head, internal resistance coefficient, and hydraulic head exponent of pump k, respectively; and ωuk is the relative speed of pump k, i.e., the ratio of working speed to the nominal speed.

Equations (3) and (5) can be unified into the general form:

Hi-Hj=qkfqk (6)

where qk is the flow rate of pipe k or pump k connecting between node i and node j; and f(qk) is defined as:

fqk=γpkqkn-1edge  k  is  a  pipe-ωuk2H0,ukqk-1+rukωuk2-mukqkmuk-1edge  k  is  a  pump (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:

qwi=kE-iqk-kE+iqk (8)

where E+(i) and E-(i) are the sets of edges (including pipes and pumps) entering and leaving node i, respectively; and qwi 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:

Awi,j=-1edge  i  leaves  node j1edge  i  enters  node j0otherwise (9)

Partitioning Aw into the following two submatrices, we can obtain:

Aw=A12 A10 (10)

where A12 is the ne×nn matrix relating the edges to the nodes with unknown hydraulic heads; and A10 is the ne×n0 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:

A11A12A12T0qH=-A10H0-qw (11)

where q is the ne×1 vector of edge flow rates; H is the nn×1 vector of hydraulic heads to be determined; H0 is the n0×1 vector of known hydraulic heads; qw is the nn×1 vector of known nodal injected flow rates; and A11 is the ne×ne diagonal matrix with the elements as:

A11k,k=f(qk) (12)

where f(qk) is the general function of pipe k or pump k.

C. Coupling Between Power System and WDN

Electricity-driven pumps are the main coupling components between the power system and the WDN. For pump k with ΔHuk and quk, the power consumption Puk is [

5]:

Puk=ρgΔHukquk/ηuk (13)

where ρ is the water density; g is the gravitational acceleration constant; and ηuk 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 QLi in (1) and (2), respectively.

III. Per-unit Modeling of Hydraulic Equations

A. Illustrative Examples of Numerical Differences

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 [0, 670]m. Other examples are pipe parameters. The pipe lengths and diameters for a typical WDN generally lie in [500, 1000]m and [100, 400]mm [

7], respectively, and a much wider range of [91, 1.56×105] for the resistance coefficients of pipes can thus be obtained with the Hazen-Williams equation.

B. Per-unitizing Hydraulic Equations

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 Table I, where the equivalent power SW is defined from the power consumption of pumps in (13).

TABLE I  Comparisons Between Variables of Power System and WDN
ScenarioParameterVariable
Power system Voltage V
Voltage drop ΔV
Current I
Impedance Z=ΔV/I
Apparent power S=VI
WDN Hydraulic head H
Head loss ΔH
Flow rate q
Resistance coefficient γ=ΔH/qn
Equivalent power SW=ρgΔHq

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 γB and SWB are determined as:

γB=HB/qBn (14)
SWB=ρgHBqB (15)

where the subscript B denotes the base value of each variable.

The resistance coefficient r depends on the individual pump, then the base value for pump k is defined as:

ruB,k=HB/qBmuk (16)

With the base values, (11) and (13) can be per-unitized, respectively, as:

A11*A12A12T0q*H*=-A10H0*-qw* (17)
Puk*=ΔHuk*quk*/ηuk (18)

where the superscript * denotes the per-unit value; and A11* is defined as:

 A11*k,k=fqk*=
   γpk*qk*n-1edge  k  is  a  pipe-ωuk2H0,uk*qk*-1+ruk*ωuk2-mukqk*muk-1edge  k  is  a  pump (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 P'uk*, which leads to:

P'uk*=Puk*SWB/SB (20)

Another point worth noting is that the hydraulic equations of a WDN should be per-unitized on a common qB for all pipes or pumps and a common HB for all nodes throughout the network.

C. Base Value Selection of Flow Rate and Hydraulic Head

Similar to the selection of SB and VB of the power system, both qB and HB can be selected quite arbitrarily. However, a well-chosen qB 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 [

7], typical ranges of pipe diameters and the economic flow rates are listed in Table II.

TABLE II  Typical Ranges of Pipe Diameters and Economical Flow Rates
Diameter (mm)Flow rate (m3/s)
100-400 0.005-0.113
>400 0.113-0.176

In the light of Table II, the following steps for selecting qB are suggested: ① calculate the average diameter of all the pipes roughly, and ② set qB to be the typical economic flow rate. For example, for a WDN with an average diameter less than 400 mm, qB is suggested to be 0.01-0.1 m3/s.

As for HB, the static heads of all the pumps shown in (5) are used by setting HB=max(H0,u1, H0,u2, , H0,um), where m is the total number of pumps.

IV. Solution to Hybrid Power-water Flow

A. N-R Method

Equations (1), (2), (11), and (13) constitute the hybrid power-water flow problem, where the power equations are expressed in per-unit value while the hydraulic equations are in actual value. This hybrid flow problem F() can be written in the general form as:

FX=0 (21)

where X=[θ  V  q  H]T 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:

Xk+1=Xk-J-1kFk (22)

where the superscripts k and k+1 denote the current iteration and the (k+1)th iteration; and J is the Jacobian matrix calculated by:

J=JeeJew0Jww=PθPVPqPHQθQVQqQH00DA1200A12T0 (23)

where Jee and Jww are the Jacobian matrices of the power system and the WDN, respectively; Jew 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:

Dk,k=nγpkqkn-1edge  k  is  a pipemukrukωuk2-muk(qk)muk-1edge  k  is  a pump (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 QLi 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 F(X) is less than the specified tolerance or the iteration count exceeds the maximum number of iterations.

B. Jacobian Matrix in Per-unit

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:

λE-J=λeEe-JeeλwEw-Jww=0 (25)

where E, Ee, and Ew are the identity matrices; and λe and λw 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 λe and therefore only λw will be discussed below. Jww can be written in (26) and λw1λw2λwn,  σ1σ2σn,  and  τ1τ2τn are the eigenvalues of Jww, B, and A, respectively, where n=ne+nn.

Jww=B+A=D000+0A12A12T0 (26)

According to the min-max theorem and interlacing theorem of eigenvalues [

8], we have:

σi+τnλwiσi+τ1λw1σ1λw2σ2σnλwni=1,2,n (27)

Considering that ① D is a positive and diagonal matrix, which means σ1=max{D(i,i)} and σn=0; and ② A(i,j)={0,1,-1} (ij), A(i,i)=0 and |A|=0 (ne>nn), 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 λwi. The 2-norm condition number of Jww, denoted by к(Jww), is proven to be max{|λwi|}/min{|λwi|}. One can decrease к(Jww) by reducing the numerical differences among D(i,i), 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:

J*=JeeJew*0Jww*=PθPVqBPqHBPHQθQVqBQqHBQH00D*A1200A12T0 (28)

where D* is the per-unit counterpart of D, which is a diagonal matrix given by:

D*k,k=nγpk*|qpk*|n-1edge  k  is  a  pipemukruk*ωuk2-mukquk*muk-1edge  k  is  a  pump (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.

C. Initialization of Unknowns

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 Vi=1 and the unknown angle θi=0. This flat start can be extended to a per-unitized WDN by setting the per-unit value of all unknown flow rates qi=1 and Hi=1. Due to the special consideration in selecting qB and HB as presented above, this initialization is usually proven to be a good start.

V. Case Studies

A. Test System Description

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.

B. Numerical Differences in Jacobian Matrix

In Table III, we compare the maximum or minimum Jacobian elements in actual and per-unit values. It can be observed that the value range of Jacobian elements is significantly reduced by per-unitizing.

TABLE III  Max/min of Jacobian Elements in Actual and Per-unit Values
ValueIEEE 33-W13IEEE 33-W27
MaxMinMaxMin
Actual value 14632.304 -21.494 576.046 -21.494
Per-unit value 21.494 -21.494 21.494 -21.494

Figure 1 shows the colormaps of the Jacobian matrix of the IEEE 33-W27 system, where the WDN is expressed in actual and per-unit values with the proposed per-unit scheme. The dimensions of Jee, Jew, and Jww are 62×62, 62×60, and 60×60, respectively, and a and b represent the coordinates of the elements in the Jacobian matrix.

Fig. 1  Colormaps of Jacobian matrix of IEEE 33-W27 system. (a) Actual value. (b) Per-unit value.

In Fig. 1(a), there is a line of bright dots in the part of Jww while the four parts of Fig. 1(b) are almost in the same color. The bright dots represent the diagonal elements of D defined in (23) and the brightness indicates that the elements of D are quite different from other elements. This is because D is determined by the hydraulic parameters of pipes and pumps while other elements of Jww are 0, 1, and -1. Moreover, since the power equations are expressed in per-unit value, the elements of both Jee and Jew are quite small compared with the elements of D in Fig. 1(a).

Both Table III and Fig. 1 verify that the proposed per-unit scheme can reduce the numerical difference among Jacobian elements significantly.

C. Condition Number of Jacobian Matrix

Table IV compares the 2-norm condition numbers of the Jacobian matrix in actual and per-unit values for the two test systems. The comparison shows that the condition numbers of the per-unit Jacobian matrix for both systems decrease significantly. Note that the condition numbers of the per-unit Jacobian matrix for the two test systems are the same. This is because the maximum and minimum singular values of the per-unit Jacobian matrix for the two systems are both determined by the IEEE 33-bus system.

TABLE IV  2-norm Condition Numbers of Jacobian Matrix
ValueIEEE 33-W13IEEE 33-W27
Actual value 4.4932×108 4.7478×106
Per-unit value 2.4641×103 2.4641×103

D. Analysis of Ill-conditioned System

For hydraulic analysis, a large variation in pipe lengths, diameters, and roughness may result in an ill-conditioned system [

6]. To investigate the impact of pipe parameters on Jacobian matrix conditioning and the proposed per-unit scheme to deal with the possible ill-conditioning, the following three cases of the IEEE 33-W27 are presented.

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.

Table V compares the 2-norm condition number of the Jacobian matrix к(J), the N-R iterations k, and the computational time T in solving the three cases. All computations are performed on a laptop, equipped with an Intel Core i5 CPU @ 1.8 GHz and 8 GB of RAM. Note that the same N-R starting point or actual tolerances for actual and per-unit values are considered. For example, the tolerance of 10-10 for hydraulic heads in per-unit value is converted to the actual tolerance by multiplying HB. Moreover, the singularity of the Jacobian matrix is also considered as the third N-R termination condition besides the tolerance and the maximum iterations mentioned above.

TABLE V  Convergence Comparison of Three Cases
CaseActual valuePer-unit value
к(J)kT (ms)к(J*)kT (ms)
1 4.7478×106 5 7 2.4641×103 5 6.5
2 5.3035×106 5 7 2.4641×103 5 7.0
3 9.8541×1014 - - 6.9513×106 6 8.1

Note:   “-” represents non-convergent.

In Table V, there are obvious differences among к(J) of the three cases, and the most obvious difference lies in case 3. Compared with к(J), the values of к(J*) are significantly reduced by the proposed per-unit scheme. Another observation is that both iteration number and computation time for the solved cases are almost the same.

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. As can be observed, the eigenvalue distribution of J* is much more concentrated than that of J. Besides, there is a very large eigenvalue of J, which results in an extremely large condition number and a convergence problem.

Fig. 2  Eigenvalue distribution. (a) J. (b) J*.

Both Table V and Fig. 2 verify the better convergence of the proposed per-unit scheme and its ability to relieve the ill-conditioning of the Jacobian matrix.

VI. Conclusion

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

1

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] 

2

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] 

3

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] 

4

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] 

5

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] 

6

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] 

7

X. Yan, “Water conveyance and distribution engineering,” in Water Supply Engineering, 4nd Ed. Beijing: China Architecture & Building Press, 1999, pp. 36-38. [Baidu Scholar] 

8

F. Zhang. (2011, May). Matrix theory. [Online]. Available: https://link.springer.com/book/10.1007/978-1-4614-1099-7. [Baidu Scholar]