Turbine modeling for steady‐state analysis in hydropower plant networks with complex layouts using a modified global gradient algorithm

Turbines are generally oversimplified in conventional algorithms for steady‐state analysis of hydropower plant (HPP) networks. Moreover, conventional algorithms are hardly suitable for serial HPPs. In this paper, a modified global gradient algorithm (MGGA) is developed by extending turbine models and optimizing turbine interpolation. First, the nonlinear model of turbines is established with the proposed variable resistance factor and added to the global matrix. Then turbine unit quantity interpolation is proposed to compute the turbine operating parameters. Finally, turbine links and pipes can be iterated synchronously, and the global matrix is solved using the Newton–Raphson method. Characteristics of MGGA are evaluated by employing two test HPPs, including rapidity of convergence, computational stability of turbine, and applicability for serial HPPs. The results show that: (a) The iterative numbers of the MGGA are 0.72 and 4.98 times on average less than that of CA1 and CA2 for the accuracy of 10−6 in 400 operating conditions, even though MGGA has an obvious hysteretic nature; (b) Only MGGA can successfully calculate the extreme cases that large head losses in HPP networks; (c) Although fixed head nodes are not present between turbine links, such as serial HPPs, MGGA still works.

some of the large hydropower bases, such as the Hindu Kush Himalayan region, due to overabundant hydropower resources, gathered ultrahigh heads, and close geological conditions, 9,10 serial HPPs with complex pipeline layouts are being inevitably developed.
The primary operating condition of HPPs is the steady state. In the design stage, the steady-state analysis of HPP networks (SSAHPPNs), one of the key features 11 of HPP operation, provides important references for (1) optimizing the hydraulic parameters of HPP networks, (2) selecting and optimizing turbines, 12 and (3) calculating the initial conditions of dynamic one-dimensional HPP simulation. 13 Hence, algorithms are employed to calculate the hydraulic parameters and the flow rate in pipes, pressures at nodes, and operating parameters of turbines under steady-state operating conditions. 14 Additionally, the operating parameters of turbines are the key parameters among all the hydraulic parameters.
Currently, HPPs with complex layouts and serial HPPs have attracted widespread concerns. However, studies on SSAHPPNs with complex HPP layouts, especially serial HPPs, are still at the initial stage. Steady-state analysis has rapidly developed for water distribution networks (WDNs), which can provide a basic reference for SSAHPPNs. Three methods are mainly employed: the loop, nodal, and pipe methods. 15 These models are usually solved using iterationbased methods, 16 such as the Hardy-Cross method, 17 the Newton-Raphson (NR) method, 18 and the linear theory (LT). 19 Note that the convergence rate of the LT could be affected by the component of pumps, valves, and so forth. 19 Currently, one of the most popular methods for simulating WDNs is the global gradient algorithm (GGA) proposed by Todini et al. 20,21 GGA takes advantage of the complete structure of the full Jacobian matrix; it can rapidly calculate the nodal heads and pipe flows in two steps and can deal with zero flow 22 and it is still being developed and improved so far. 21,23,24 Unfortunately, it has not been utilized for turbine models. A High-order global algorithm for pressure-driven modeling of WDNs is proposed based on increasing the order of convergence. 25 The modelizations of various valves in different methods have been investigated in many studies. [26][27][28][29] Nevertheless, the turbine model for steady-state analysis has not been investigated in detail. Furthermore, the strong nonlinear characteristics of turbines are complex in SSAHPPN modeling . [30][31][32] Quite a number of meaningful investigations have been conducted for turbine modeling in the past decades. An ANN method of the correlations for efficiency using operating parameters is developed by Kumar et al., 33 which is applied to evaluating the performance of turbines. An Adaptive-Neuro-Fuzzy Inference system is proposed by Weldcherkos et al., 34 which can help to speed up the performance of the automatic generation control. Flow losses in a turbine are modeled with the entropy production theory by Wang et al., 35 which can be feasible to explore where the energy dissipation happens in turbine passage under different operating conditions. Although many meaningful works have been carried out for turbine modeling. However, only a few relevant studies on turbine modeling for SSAHPPPNs have been reported. Zhao et al. calculated the initial variables of a real-time accurate equivalent circuit model using a simplified algorithm for tree networks. 4 Yang et al. introduced TOPSYS, a transient simulation software, 36,37 including a calculation module of the steady-state analysis using a simplified fixed point iteration algorithm which has been applied to many studies. 37,38 However, the two abovementioned conventional algorithms are only applicable to single HPP networks, and they are simply reviewed herein. The turbine model is oversimplified, then SSAHPPNs are computed without a complete solve system in these two conventional algorithms. Moreover, to the best of our knowledge, an algorithm of steady-state analysis for serial HPPs has not yet been reported. The above literature review and analysis indicate that the features of the existing works and growing trend are as follows for detail: 1. With the increasingly complex layouts of HPPs, especially serial HPPs, 39,40 an algorithm needs to be developed including the turbine model for the refined SSAHPPNs with complex HPP layouts. 2. Turbines are oversimplified, and the system is complexly solved by partitioning pipes and turbines in conventional SSAHPPN algorithms. Moreover, conventional algorithms are only applicable for the topological layouts of the case always existing fixed head nodes between turbine links, and they are hardly suitable for serial HPPs with complex turbine layouts. 3. GGA is the most widely used and extended algorithm among all the algorithms. To the best of our knowledge, a refined model of turbines for steady-state analysis, the key component in SSAHPPNs, has been rarely reported in the literature. GGA could be extended and studied in depth.
Aiming at these features, based on the traditional GGA, an improved algorithm, referred to as the modified GGA (MGGA), is proposed for SSAHPPNs herein. The nonlinear model of turbines is investigated and included in the GGA global matrix. Hence, the complete system can be solved using the NR method. The characteristics and advantages of MGGA are summarized via case studies. MGGA can be innovatively applied to the steadystate analysis of serial HPPs. Moreover, for the cases of single HPP, compared with conventional algorithms, it affords great computational stability in extreme cases and a fast convergence rate in a majority of operating conditions. The overall process of this study is shown in Figure 1.
The remainder of this paper is organized as follows. In Section 2, based on the fundamentals of HPP network hydraulics, the turbine model and variable resistance factor are proposed, and then, the MGGA derivation is presented. Moreover, the conventional algorithms are simply mentioned. In Section 3, two test HPPs are employed to investigate MGGA properties. Then, the characteristics of MGGA are summarized. In Section 4, MGGA is compared comprehensively with conventional algorithms, and the limitations of MGGA are listed. In Section 5, the conclusions and recommendations are presented.

| Fundamentals of HPP network hydraulics
Generally, an HPP network is represented by links that convey water between nodes that connect links to one another. Figure 2 displays a complex HPP network, describing the essential elements in HPP systems. Short forms are used to describe serial numbers of pipes and turbines according to the topology of an HPP network: for example, pipe 1, Node 2, and turbine 5 are called P1, N2, and T5 in the HPP network of Figure 2, respectively.
Serial HPPs are special cascade HPPs that can be explained as two or more HPPs that are hydraulically connected end-to-end by small water reservoirs, called regulating reservoirs. Owing to the close geographical and hydraulic connection, the operation of serial HPPs needs to be coordinated. In summary, it is a novel form of HPPs with a complex layout of turbines and pipelines. Figure 3 presents a simple diagram of serial HPPs.

| Definitions and notations
Consider an HPP network comprising p pipe links, t turbine links, n unknown head nodes, and n 0 fixed-head nodes. To easily describe the topology of an HPP network, turbines are simplified and considered as particular links herein, instead of as compositive energy elements in conventional HPP simulation software TOPSYS and SIMSEN. Subsequently, the HPP network comprises P (p + t) links and N (n + n 0 ) nodes.
The kth pipe of the network is characterized by its diameter D k , length L k , and resistance factor r k . The lth F I G U R E 1 Overall process of the current study MA ET AL.
| 1253 turbine of the network is characterized by its runner diameter D 1 , characteristic curves, GVO α 0 , and rated speed n r . Additionally, the ith node of the network has two properties: its nodal demand q i and its elevation head z i . head nodes and fix-head nodes, A 12 is partitioned as follows: 12 12 10 (2) where A 12 is a [P; n] matrix relating the links to the unknown head nodes, and A 10 is a [p; n 0 ] matrix relating the pipes to the fixed head nodes. For clarity, the transposes of matrices are expressed simplified: e.g. A 21 represents A T 12 and A 01 represents A T 10 herein.

| Basic equations of pipe networks
The energy and continuity equations are the basic equations 19 employed to describe the steady-state flows and heads in the HPP networks. The energy equation describes the head loss through each pipe or head for each turbine: where subscripts i and j represent the start and end nodes of link k, and g(Q k ) represents the head loss due to resistance at pipe k or head for energy conversion at turbine k as a function of its flow. The turbine head is introduced in Section 2.1.4. The head loss through a pipe is where r k represents the resistance factor, characterizing the frictional and local head losses because the local head losses are the major losses in HPP networks that cannot be omitted. Additionally, α = 2 while using the Darcy-Weisbach equation in HPP networks because it is convenient to combine the frictional and local head losses.
In addition to Equation (3), the mass continuity equation describes the sum of inflows and outflows as zero at each unknown head node i: where Q ki,j represents the flow in the link k i,j connected node i; n i represents the number of links connected to node i; and q i represents the nodal demands, which is generally zero in HPP networks.
In an HPP network with P links and N (n + n 0 ) nodes, Equations (3) and (5) can be expressed as a global matrix, which is adopted by GGA: Here, A 11 represents a diagonal matrix whose elements are defined as for k [1, P].

| Nonlinear model of turbine
Compared with that in a WDN, the operating parameters of turbines are more crucial than the flow of pipes and the head of nodes in an HPP network. Hence, the refined modeling of turbines needs to be investigated. The most remarkable characteristics of turbines are their strong nonlinear characteristics. The discharge of turbine Q t comprises the nonlinear functions of head H t , unit speed n 11 , and GVO α 0 , which is described by the characteristic curves without an accurate function expression. 30 The unit quantities in characteristic curves are defined as follows: where Q 11 and n 11 represent the unit discharge and unit speed, respectively. Characteristic curves obtained via model experiments are generally sparse, and incompletely measured curves from the manufacturer have an adverse effect on the steady-state analysis. Furthermore, unqualified data points in measured curves need to be corrected to decrease the measuring error and improve the model accuracy. Hence, characteristic curves need to be preprocessed via internal interpolation. For Francis turbines, preprocessing can be performed using a BP method. For pump turbines, it can be performed via an improved Suter-BP method. 4 After the above preprocessing, the introduced refined characteristic curves can be used to solve the nonlinear function of turbine models and calculate the variable resistance factor r t .

| Derivation of the modified global gradient algorithm including turbine
Owing to the strong nonlinear characteristics of turbines, turbine links are generally disconnected from the function matrices in conventional algorithms. Consequently, the topological structures of HPP networks are incomplete in the function matrix of conventional algorithms. To complete the topological structures of HPP networks and reflect the nonlinear characteristics of turbines in the function matrices of SSAHPPNs, MGGA, which includes turbines, is developed herein.
To highlight the turbine links in the global matrix, some of the matrices are modified as follows: n] matrix relating the pipes to the unknown head nodes, and A 13 is a [t; n] matrix relating the turbines to the unknown head nodes. Substituting these into Equation (6), we get 2.2.1 | Variable resistance factor of turbine Q t , H t , and r t of a turbine link vary with the different hydraulic conditions. This is because when the water level of upstream and downstream reservoirs changes or the demand of the power grid changes, the power output and discharge of turbine Q t should be changed by adjusting GVO. To reflect the inherent flow characteristics of turbines, avoid the repeated conversion of unit quantities and parameters, and improve the interpolation calculation efficiency, the function of r t and unit quantity in the characteristic curves need to be preferably found.
When a turbine operates in the grid-connected mode, the rotational speed of the turbine is the rated speed n r . Also, the kinetic energy correction H Δ should be considered: where β t denotes the kinetic energy correction factor of the turbine, which is a fixed constant; α p and α s denote the kinetic energy correction factor of the spiral case outlet and draft tube inlet, respectively; A p and A s denote the kinetic energy correction factor of the spiral case outlet and draft tube inlet, respectively. Then, r t can be expressed by the unit quantities and diameter D 1 based on Equation (8). The variable resistance factor r t of a turbine link is proposed as follows: The crucial process involves solving the nonlinear model of the turbine; hence, the variable resistance factor r t can be calculated using Equation (14). Q t can be calculated from the refined characteristic curves, D 1 , α 0 , H t , and n r . Herein, the solve process is called the turbine unit quantity interpolation (TUQI).
For the special case of an outage, the GVO of the turbine is zero, that is, no flow exists in the pipeline. To simplify the calculation and retain the original network topology, a huge resistance factor (r t = 10 5 ) is employed herein to replace the zero flow model of the outage turbine. Notably, H t represents the pressure difference between the front and rear nodes, instead of the turbine head for the outage case.

| MGGA and its NR iteration
The following equations for the flow and head corrections at iteration τ of the process are afforded when directly applying the NR method to Equation (12): Where D pp and D tt represent the Jacobian of A pp Q p and D tt Q t , respectively, and k represents the kth diagonal element.  Moreover, dQ p and dQ t can be solved based on Then, dH can be solved by Equation (23). The iterative formulation of MGGA, Equations (24)-(26), can be found by substituting f f f , , and from Equations (19)-(21).
where I pp and I tt are the identity matrices of size [p; p] and [t; t], respectively. The final iterative formulation of MGGA can be produced by eliminating the identity matrices and transposing, as shown in Equations (28)- (30). Figure 4 shows the iterative steps of MGGA.

| Conventional algorithms of SSAHPPNs for comparison with MGGA
Two conventional algorithms, introduced in Ref. 4 For each turbine, the "upstream reservoir-turbinedownstream reservoir" pipeline is divided from the HPP network. Therefore, it can be simplified as a tree network. It is the simplest algorithm for simulating SSAHPPNs based on the special network layout of HPPs. The iterative steps of CA1 are as follows: Step 1: Initial H t (0) is selected as the head difference between the upstream and downstream reservoirs: where l represents the lth turbine.
Step 2: Calculate Q t τ ( ) using the characteristic curves; Step 3: Meet the mass continuity conditions: Step 4: Update H t τ ( +1) using the energy conditions: Step 5: Evaluate whether the iteration meets the accuracy target. If the iterative accuracy is satisfied, go to Step 6; otherwise, return to Step 2.
Step 6: Final H can be computed as

| CA2: Simplified fixed point iteration algorithm
A function matrix of CA2 is created using the mass continuity conditions. The matrix equations are solved by the NR method. Imperfectly, for turbine links, H t is not considered in the function matrix. Hence, the topological structure of HPP networks is incomplete in the function matrix. CA2 is adopted in the transient simulation software TOPSYS. The iterative steps of the CA2 are as follows: Step 1: Initial H t (0) is selected, similar to that in CA1.
The initial H (0) can be selected as a unit vector; Step 2: Calculate Q t τ ( ) using the characteristic curves.
The initial Q (0) is selected similar to Q t τ ( ) ; Step 3: Update H (τ+1) and H t τ ( +1) : where Step 4: Update H t τ ( +1) using the characteristic curves; update Q (τ+1) as follows: Step 5: Evaluate whether the iteration meets the accuracy target. If the iterative accuracy is satisfied, end the iteration; otherwise, return to Step 3.

| Test HPPs
The MGGA and conventional algorithms are compared using two typical HPP networks, including a single and a serial HPP network, which are shown in Figures 5 and 6, respectively.
Notably, two cascade HPPs are selected as serial HPPs and connected by a hypothetical regulating reservoir. Moreover, some parameters of turbines in the JP 2 HPP are modified to satisfy the continuity equations and comprehensively consider the operating efficiency of serial HPPs. The crucial parameters of HPPs and the modifications are listed in Table 1. For node characteristics, nodal demands q of two HPP networks are zero vectors. Tables A2 and A3 list the characteristic parameters of pipes. The characteristic curves of turbines are shown in Figure 7.

| Calculations
Three tests were carried out using test HPP networks. The first one is set to assess the convergence rate of MGGA in different required accuracies, also the correctness of MGGA is verified; The second one is intended to illustrate the computational stability of the turbine, especially in extreme operating conditions; The last one is designed to verify the applicability of MGGA for the serial HPP networks. The executions have been performed on a desktop computer with an Intel Core i7-9700 CPU at 3.00 GHz and 16.0 GB RAM, running Windows 10. All the codes are run on the MATLAB platform R2021a and have the same optimization options in all cases.
Each algorithm is commenced with the same initial set of conditions: the initial discharge of turbines is set as the rated discharge, the flows of other pipes are set as 1 m 3 /s, and the initial heads of nodes are set as 1 m.
Accuracy is defined as the maximum relative deviation between the computed flow in iteration τ and iteration τ − 1 in 1 s −1 . The relative error is defined as absolute errors relative to the average calculations of all the calculated algorithms. The absolute errors are defined as the absolute value of the difference between the corresponding algorithms and average values. Relative value is defined as the iterative number or run time of other algorithms relative to MGGA.

| Rapidity of convergence
Using the LDL HPP as an example, various GVO combinations of turbines are considered to comprehensively demonstrate the iteration and convergence properties of MGGA, also the correctness of MGGA can be verified. GVO commences at 0.05 and increases in steps of 0.05 until it reaches 1. Since two turbines are present in LDL HPP, a total of 400 operating conditions are calculated. The water levels of the upstream and downstream reservoirs in each operating condition are NPL. To compare both low and high accuracy requirements, the required accuracies are set as 10 −2 , 10 −4 , and 10 −6 .  Moreover, the reference results are provided to investigate the influence of turbine interpolation. The correct r t is directly assigned to the turbine links to calculate the reference results, the turbine interpolation step is omitted, and the reference is solved by GGA.
The global iterative numbers and run time of each algorithm are shown in Figure 9; moreover, the errors after the first iteration are provided to investigate the convergence properties, as shown in Figure 10. Clearly, for reference, without the turbine iteration steps, the fastest convergence is afforded. The iterative numbers of MGGA are slightly less than those of CA2 and much less than those of CA1, which can see from the average and median values of iterative numbers in Figure 9A. For the run times, MGGA runs the fastest, followed by CA2 and CA1, as shown in Figure 9B.
The detailed iteration features of different algorithms are listed in Table A4. According to the calculation of the above 400 operating conditions, MGGA gets the best iterative convergence performance among the three algorithms; Moreover, with increasing required accuracy, the rate advantage of the MGGA becomes more significant. For the iterative numbers: Compared with the reference results, the iterative numbers of the MGGA increase by 0.77, 1.03, and 1.20 times on average for the three required accuracies after the addition of the turbine interpolation step, which is also 0, 0.29, and 0.72 times on average less  than that of CA2 and 1.07, 2.96, and 4.98 times on average less than that of CA1. For the run time: Compared with the reference results, the run time of the MGGA increased by 2.41, 3.35, and 5.18 times on average for the three required accuracies after the addition of the turbine interpolation step, which is also 0.59, 1.94, and 1.94 times on average less than that of CA2 and 2.94, 7.00, and 9.06 times on average less than that of CA1.
The iteration and convergence properties are described in detail in Section 3.3.

| Computational stability of turbine
The extreme case means a large head loss is modeled in the HPP networks in this research: such as when a gate is almost shut in the pipeline or shaft, or a ball valve is operated with a small GVO or the error modeling of HPP networks by technicist. It may present a challenge for the computational stability of turbines in algorithms.
Two extreme cases are calculated in LDL HPP networks to illustrate the computational stability of turbines. Only resistance factors are changed: r of P12 and P14 increase to 10 −3 in OC5 and OC6, respectively, which is a challenge to turbine interpolation of T5 in OC5 and T12 in OC6. Other parameters of the two cases are the same: water levels of the upstream and downstream reservoirs are still in the NPLs, both the GVOs of T5 and T12 are set as 0.9153, and the required accuracy is set as 10 −6 .
The calculations are shown in Table 2. Only MGGA can successfully perform the computation, and both CA1 and CA2 fail. The calculations of MGGA can be deemed correct according to the comparison in Section 3.2.1.
Combinations of different initial conditions are calculated to illustrate that the failure of CA1 and CA2 is independent of the initial conditions. Six representative initial conditions are selected for exposition: Both the Q of T5 in OC5 and that of T12 in OC6 are set as 1, 101, 202, 303, 404, 506 m 3 /s respectively; moreover, the H of all nodes are set as the correct values in CA1 and CA2 to reduce the difficulty of iteration; and other parameters remain unchanged. The iterative processes of H of the turbine link are shown in Figure 11.
The crucial process in the steady-state simulation for HPP networks is turbine interpolation. For conventional algorithms, turbine links are generally disconnected from function matrices, and node heads are calculated using function matrices. Subsequently, the turbine heads are calculated. Finally, the discharges of the turbine are calculated by interpolation using characteristic curves. Mostly, the parameters and layout of the HPP networks satisfy the conventional turbine interpolation condition, and the iterative process continues. However, in extreme cases OC5 and OC6, the conventional turbine interpolation is observed to be almost unstable because the turbine head during the iterative process is calculated to a negative value. Consequently, the computation cannot be performed using Equation (8). CA1 and CA2 should be improved to meet the challenge of extreme cases.
MGGA overcomes the abovementioned disadvantage: the turbine links are described using the unit quantity and variable resistance factor, and hence, the topological structures of HPP networks in the global matrices are completed. Moreover, TUQI avoids the transition between the actual and unit parameters of the turbine, which improves the computational stability of the turbine.

| Applicability for serial HPPs
For serial HPPs, at least two serial turbine links are layout in a pipeline, no fixed head nodes are present between turbine links additionally. Hence it is a great challenge for SSAHPPNs. As shown in Section 3.2.2, CA1 and CA2 are not Applicable to the cases of large head losses that existed in the HPP networks; Predictably, they are also not applicable to the serial HPPs. Similarly, the applicability of MGGA for serial HPPs should be re-verified.
To verify the MGGA calculations and without a reported efficient algorithm for serial HPP networks, CA2 is simply improved and nested inside the dichotomy, which is called "CA2-D" herein. The required accuracy of the dichotomy is set as 10 −2 . First, the water level of the regulating reservoir is initialized as the  average value of the water level of the upstream and downstream reservoirs. Then, the steady states of JP 1 and JP 2 are separately simulated using CA2. Based on the continuity equations at the regulating reservoir, the water levels of the regulating reservoir are updated via dichotomy. Subsequently, repeated iterations are performed. Finally, if the continuity equations are satisfied, the iteration is completed. Three typical operating conditions are considered to simulate the steady state of serial HPPs: outage, symmetric operation, and asymmetric operation (Table 3). Operating condition 1 and operating condition 2 are corresponding to a high and a medium power output of symmetric operations, respectively; operating condition 3 and operating condition 4 are corresponding to asymmetric operations; for operating condition 4, turbines T19 and T24 are in outage additionally. The water levels of the upstream and downstream reservoirs in each operating condition are at normal poor levels (NPLs), and the required accuracy is set as 10 −6 . Figure 12 displays the computed head and discharge of each turbine for the four operating conditions calculated by MGGA and CA2-D. Tables 4 and 5 display relative errors and iterative parameters for the four operating conditions calculated by MGGA and CA2-D algorithms, respectively.
1. For result accuracy, the Results of CA2-D and MGGA are basically the same, which shows that MGGA is reliable. The maximum relative errors of discharge and head are 6.26% and 0.36%, respectively. The maximum relative errors of flow are corresponding to the turbine in the outage (T4), which are equivalent to the absolute errors of 0.05 m 3 /s. It gets a large relative error because the calculated flow is almost zero. Eliminate turbines in the outage, the maximum relative error of discharge is 0.53%. The errors might result from the different GVO interpolation of characteristic curves in the two algorithms; Furthermore, the required accuracy at the regulating reservoir (N) of CA2-D is more imprecise than that of MGGA. 2. Note that CA2 cannot be directly used in the steadystate analysis of serial HPPs. However, using preliminary algorithm nesting, the improved CA2-D can be utilized for the steady-state analysis of serial HPPs.
The results show that CA2-D affords poor  In summary, the steady state of serial HPPs can be simulated using MGGA, but it is difficult using conventional algorithms. Moreover, MGGA affords a fast convergence property, and its convergence rate is little affected by the model of the outage turbine.

| Explanations and summaries
1. The iteration and convergence are clarified from three aspects: the principle of the algorithm, the required accuracy, and the error after the first iteration: • Compared with the reference results, for the different required accuracies, the iterative numbers of MGGA increase. This is because the convergence behavior of MGGA is affected by the nonlinear characteristic curves of the turbine after the addition of the turbine interpolation step. • MGGA shows an obvious hysteresis nature in the first iteration. In 68% of the operating conditions, the error of the first iteration is larger than that of the conventional algorithms, and the error after the first iteration increases with decreasing GVO. After considering the turbine links in the global matrix, the turbine links and pipes are synchronously iterated. For MGGA, since the internal strong convergence of the characteristic curve is negatively affected by the iteration of other pipes, the error after the first iteration is large. For CA1 and CA2, the error of the first iteration is small due to the separate calculation of the turbine interpolation. • Generally, MGGA affords the best convergence rate, inheriting the characteristic of fast convergence of GGA. As the required accuracy increases, the rate advantage of MGGA becomes more obvious under the premise of the hysteretic nature. • In summary, the iterative numbers of MGGA are affected by the required accuracy and convergence rate. When the required accuracy is low, the calculations are mainly affected by the hysteretic nature, increasing the iterative numbers. As the required accuracy increases, the effect of the convergence rate increases, and the rate advantage of MGGA becomes more significant. 2. For the computational stability of the turbines, the laws of turbine interpolation for MGGA and conventional algorithms are illustrated, and the extreme cases are analyzed.
• The results show that only MGGA can successfully perform the computations, and both CA1 and CA2 fail unrelated to the set of initial conditions, indicating that the possible error of turbine interpolation is avoided using MGGA. The topological structures of HPP networks in global matrices are completed by describing the turbine link as the unit quantity and variable resistance factor, and the turbine links and pipes are synchronously calculated. Hence, MGGA affords the best computational stability of the turbine among the three algorithms investigated. • The high computational stability of the turbines increases the error tolerance rate of the steady-state simulation and avoids the error in the computation process, which makes it easier for technicians to examine the modeling of HPP networks and simulate the steady state. Furthermore, more extreme conditions of HPP networks can be calculated using MGGA. 1. For conventional algorithms, the convergence rate is mainly affected by turbine interpolation. They are the "internal iteration of function matrix with external turbine interpolation" algorithms with incomplete topological structure in the function matrices. 2. For MGGA, the turbine links are described as unit quantities and variable resistance factors and added to the global matrix. Additionally, TUQI is used for turbine interpolation. Therefore, the abovementioned disadvantages are resolved, and MGGA affords great computational stability and a fast convergence rate. Moreover, based on the complete topological structure in the global matrix, MGGA can be used for a wide application scope, such as serial HPP networks. The comparison of MGGA and conventional algorithms is shown in Figure 13 and Table 6.

| Limitations and future work
This study has several limitations, and could be extended in the following aspects: 1. The modeling of outage turbines can be further developed. The outage turbine might be removed from the global matrix without disconnecting the HPP networks by a method of detecting and labeling, for reducing the matrix size and accelerating the solving speed.
F I G U R E 13 Essential comparison of MGGA and conventional algorithms. MGGA, modified global gradient algorithm.
2. Further work can be conducted on more hybrid networks comprising open channels and pipes. Open channels may be present in the bridging links of serial HPPs and sloping ceiling tailrace tunnels; Also more regulating reservoir types and Pelton turbines could be investigated in-depth. Hence, the application scope of MGGA could be widened. 3. MGGA could be extended to developing unsteady flow algorithms. For the pipe modeling, the spatial discretization of pipes is not required both for MGGA and rigid water column models. Hence, rigid water column models and refined turbine modeling could be investigated and improved based on MGGA, for simulating slow transients of HPP networks effectively in the future.

| CONCLUSIONS AND RECOMMENDATIONS
In this study, MGGA, an algorithm that can be employed for the steady-state simulation of HPPs with complex layouts, especially serial HPPs, is developed. First, based on the traditional GGA, the nonlinear model of turbines is investigated, which is described as unit quantities and variable resistance factors. Then, turbine links are added to the global matrix, and the process of turbine interpolation is optimized using TUQI. Two test HPPs are employed to evaluate MGGA, including the rapidity of convergence, the computational stability of the turbine, and the applicability of serial HPPs. Compared to conventional algorithms, the topological structure of HPP networks in the matrix is complete. Thus, MGGA can be used for a wide application scope, and it affords great computational stability and a fast convergence rate.
The key findings and contributions are as follows: 1. The major contribution of this study is the modeling of the turbine link and the improvement of the topological structure of HPP networks in the global matrix: • MGGA affords a fast convergence rate: The process of turbine interpolation is optimized using TUQI, avoiding the repeated transition between the actual and unit parameters of the turbines. Then, with the complete topological structure of HPP networks in the global matrix. • MGGA affords great computational stability, which can be used for simulating extreme operating conditions. The turbine links are described as variable resistance factors; Then, turbine links are added to the global matrix. The turbine links and pipes are synchronously iterated. • Wide application scope of MGGA is elucidated for not only HPPs with a complex layout but also serial HPPs. Although fixed head nodes are not present between turbine links, such as in the layout of "turbine-regulating reservoir -turbines," MGGA can still simulate SSAHPPNs. 2. For engineering applications, MGGA can be employed for the steady-state simulation of HPPs with complex layouts, especially serial HPPs.
• MGGA retains the characteristic of fast convergence of GGA. Compared with conventional algorithms, the improvement of the convergence rate saves considerable time during the hydraulic optimization of HPPs and the simulation of hydraulic transients. • Moreover, extreme operating conditions can be calculated, making it easier for technicians to examine the modeling of HPP networks and simulate the steady state, and the calculation efficiency is effectively improved.
In addition, this study could be extended in these aspects: refined modeling of an outage turbine, modeling of channel and Pelton turbines, and simulation of slow transients.