## ABSTRACT

In a complex pressurized water diversion project (WDP), the combined optimal operation of multiple hydraulic facilities is computationally expensive owing to the requirement of massive mathematical simulation model runs. A parallel multi-objective optimization based on adaptive surrogate model (PMO-ASMO) was proposed in this study to alleviate the computational burden while maintaining its effectiveness. At the simulation level, an adaptive surrogate model was established, while a parallel non-dominated sorting genetic algorithm II (PNSGA-II) was utilized at the optimization level. Taking the successive shutdown of pumps as the operating process, the PMO-ASMO was applied to a complex pressurized diversion section of the Jiaodong WDP in China, and the results were compared with those obtained by NSGA-II and PNSGA-II. The results showed that the time consumption of PMO-ASMO was only 9.97% of that acquired by NSGA-II, which was comparable to that of PNSGA-II, in the case of 10-core parallelism. Moreover, compared with PNSGA-II, PMO-ASMO could find the optimal and stable Pareto front with the same number of simulation model runs, or even fewer runs. These results validated the effectiveness and efficiency of the PMO-ASMO. Therefore, the proposed framework based on multi-objective optimization is efficient for combined optimal operation of multiple hydraulic facilities.

## HIGHLIGHTS

A combined optimal operation of multiple hydraulic facilities within complex pressurized water diversion projects is achieved.

The combined optimal operation model of multiple hydraulic facilities coupled with the simulation model is established.

The algorithm that combines the adaptive surrogate model and parallel computing is provided.

The effectiveness and efficiency of the proposed method are confirmed by an engineering case.

## INTRODUCTION

Water diversion projects (WDPs) have emerged as effective measures to address the uneven distributions of water resources, thereby improving the ecological environment and promoting regional economic development (Ma *et al.* 2020). With the worldwide development of WDPs, the importance of their operational strategies is constantly growing (Zhao *et al.* 2019). In particular, for complex pressurized WDPs with long water pipelines and multiple hydraulic buildings and facilities, the optimal combined operation of multiple hydraulic facilities is essential for ensuring a safe transient process under highly nonlinear hydraulic-mechanical coupling characteristics.

An efficient multi-objective optimization algorithm (MOA) is an effective approach to addressing the problem (Zhang *et al.* 2017). However, compared with the regulation of a single hydraulic facility of a simple WDP (Skulovich *et al.* 2016; Xu *et al.* 2018; Lai *et al.* 2019; Zaki *et al.* 2019; Rezghi *et al.* 2020), the transient processes under the combined regulation of multiple hydraulic facilities in a complex pressurized WDP are intricate (Cui & Guo 2023). The coupling relationships among multiple hydraulic facilities have an important effect on water hammer wave. How to determine the coupling relationships based on the certain safe operation requirements? Moreover, since the multi-objective optimization (MO) problem requires extensive simulations to obtain the global Pareto optimal solutions (Gong *et al.* 2016a), and the mathematical simulation model of complex pressurized WDPs requires tens of minutes or even hours to run once, it is difficult to efficiently obtain a multi-objective Pareto frontier for the combined optimal operation of multiple hydraulic facilities. How to reduce the computational cost? The above issues are the obstacles for the combined optimal operation of multiple hydraulic facilities in complex pressurized WDPs.

For the combined operation of multiple hydraulic facilities, Lai *et al.* (2020) employed the combined optimization of guide vane closing laws for two units under successive load rejection conditions in a pumped storage hydropower station. Zhao *et al.* (2019) investigated the optimal coordinated operation of the ball valve and guide vane based on a reference vector-guided evolutionary algorithm. Hou *et al.* (2019) studied the effects of time interval on startup quality during the successive startup of pumped storage hydropower stations by a multi-objective particle swarm optimization algorithm (MOPSO). Wan & Li (2016) proposed a joint pump-valve regulation scheme for a series pump-valve system, considering the effect of the time interval between the valve and pump during the startup process. Yu *et al.* (2015) proposed a unique closure law for the joint regulation of main inlet valves and wicket gates under the load rejection process for a pumped storage hydropower station. Zhang *et al.* (2022) adopted an improved multi-objective gravitational search algorithm (MOGSA) to investigate the joint optimal control of the valve and guide vane under the load rejection condition of hydroelectric units. Lei *et al.* (2021) optimized the startup strategy for a hydroelectric system based on MOPSO. Cui & Guo (2023) studied the multi-objective control of transient process of hydropower plant with two turbines sharing one penstock under combined operating conditions. The most favorable superposition time and the most unfavorable superposition time under COCs were determined. However, few studies have been conducted on the combined optimal operation of multiple hydraulic facilities for complex pressurized WDPs with long water pipelines and multiple hydraulic buildings and facilities.

Parallel computing and surrogate-based optimization are efficient measures for reducing computational time. Parallel computing is achieved by dividing a complex father task into multiple subtasks that can be solved parallelly by different processors. With the rapid advancement of electronic technology, the majority of computers are equipped with multicore processors and provide the necessary platform for realizing parallel computing. Additionally, numerous companies have succeeded in developing simple yet powerful parallel frameworks, thereby demonstrating essential advanced software platforms and conditions (Feng *et al.* 2018a). In MO problem solving, parallel processing technology was combined with optimization algorithms, including multi-objective genetic, particle swarm, and shuffled frog leaping algorithms, that efficiently address the optimization problem in reservoir operations (Sun *et al.* 2016; Feng *et al.* 2018b; Yang *et al.* 2022) and dispatching problem in battery storage systems (Grisales-Noreña *et al.* 2020). Surrogate models, also known as metamodeling, can accurately mimic complex models. Hence, surrogate-based optimization can effectively reduce the computational burden by replacing previous simulation models with cost-effective statistical surrogates. It has been widely applied in the field of groundwater optimization problems, the design and optimization problems of water distribution systems, and so on (Johnson & Rogers 2000; Behzadian *et al.* 2009). Recently, MO methods based on adaptive surrogate model have emerged for the optimization of the design and operation of WDPs. For example, Bazargan-Lari *et al.* (2013) adopted the MO model and Bayesian networks to optimize the closing rule for the valve of a simple reservoir-pipe-valve system, thereby avoiding the time-consuming simulation and optimization processes. Tong *et al.* (2020) built three surrogate models by analyzing the highly nonlinear relationship between the external characteristic values and essential design variables of centrifugal pumps. The three surrogate models were separately used to optimize the essential design variables in combination with the non-dominated sorting genetic algorithm II (NSGA-II). Gong *et al.* (2016a) proposed an MOA based on adaptive surrogate model, which was applied to parameter estimation and reservoir optimization operations in large and complex geophysical models. Wu *et al.* (2023) built a surrogate model using Gaussian process regression (GPR) and used NSGA-II to optimize the impeller and diffuser. With regard to the combined optimal operation of multiple hydraulic facilities, few works have been conducted using optimization algorithms combined with an adaptive surrogate model. Additionally, there exist few reports on the integration of parallel computing and adaptive surrogate models.

This study aims to utilize a novel approach called PMO-ASMO, which integrates parallel computing, MO, and an adaptive surrogate model, to obtain the combined optimal operation schemes of multiple hydraulic facilities in complex pressurized WDPs. The novelties and innovations of the study can be concluded as follows: (1) A MO model for the combined operation of multiple hydraulic facilities in complex pressurized WDPs was established with full consideration of the coupling relationship between hydraulic facilities. (2) An efficient method for the optimization model based on the parallel computing, adaptive surrogate model and optimization algorithm was proposed.

In Section 2, the mathematical simulation model for complex pressurized WDPs was established, and transient interference analysis of coupling relationships between multiple hydraulic facilities was analyzed, then a MO model for the combined operation of multiple hydraulic facilities was built, and the efficient method called PMO-ASMO was proposed. In Section 3, the Pareto frontier of the combined optimal operation of multiple hydraulic facilities achieved by PMO-ASMO was compared with those obtained by the classical evolutionary algorithm to verify the effectiveness and efficiency, taking a complex pressurized diversion section of Jiaodong Water Diversion Project in China as example. In Section 4, the conclusions were given.

## METHODOLOGY

(1) Mathematical simulation model for complex pressurized WDPs: The general layout of complex pressurized WDPs comprises a long main pipe with multiple series-connected hydraulic buildings and facilities, including pumping stations, high-level water tanks, branch pipes, surge tanks, valves, and regulating tanks. Further, the pumping station generally comprises several pump-valve units. Each unit is formed of a valve and a pump in series, which is joint to the main pipe via a branch pipe. Hence, to establish a mathematical simulation model, the transient flow in the pipes is described using governing equations, while hydraulic buildings and facilities are treated as boundary conditions that are coupled to the governing equations. The method of characteristics (MOC), which is a simple, numerically efficient, and the most used method, was adopted to solve the equations.

(2) Transient interference analysis of coupling relationships between multiple hydraulic facilities: Given the general layout of complex pressurized WDPs, the successive shutdown of pumps was adopted as the operating process in this study. The mutual transient interference of the coupling relationship among multiple hydraulic facilities is explored using a mathematical simulation model to identify the adjustable parameters (inputs

*P*) and response functions (outputs*O*) of the model when facing an actual problem.(3) Multi-objective optimization: When adopting the PMO-ASMO, firstly, multiple GPR surrogate models are constructed for different response functions. To establish the surrogate model, sampling techniques are utilized, and their initial sampling is arranged using the Good Lattice Points (GLP) method with ranked Gram-Schmidt (RGS) decorrelation. A multicore parallel technique using the Fork/Join framework is integrated into the initial sampling process. Secondly, the PNSGA-II algorithm, which incorporates parallel technology and the NSGA-II algorithm, is employed in the GPR to optimize the surrogate models and obtain Pareto optimal points. When referring to NSGA-II and PNSGA-II, the NSGA-II is utilized in a serial and parallel manner to drive the mathematical simulation model of the complex pressurized WDPs to obtain Pareto optimal points, respectively. Finally, to evaluate the validity and advantages, the PMO-ASMO was compared with the PNSGA-II and NSGA-II.

(4) Pareto frontier: The results of the aforementioned algorithms were analyzed to identify appropriate combined operating schemes for multiple hydraulic facilities.

### Mathematical simulation model for complex pressurized WDPs

A specific mathematical description of each component of the WDPs is the foundation for optimization. In this section, a mathematical simulation model that includes the pipe, pump, and valve is presented.

#### Pipe

*V*≪

*a*, therefore Equations (1) and (2) can be simplified to Equations (3) and (4), which are further converted into the following difference equation (Equations (5) and (6)) by the MOC (Wylie

*et al.*1993; Wang

*et al.*2023):where

*H*is the instantaneous piezometric head,

*Q*denotes the instantaneous discharge,

*V*represents the average speed,

*i*is the section number along the pipeline,

*j*stands for the time in the algebraic water hammer,

*α*is the pipe slope,

*g*is the gravitational acceleration,

*a*is the wave speed, and

*D*and

*f*are the diameter and friction factor of the pipe, respectively.

*C*,

_{P}*C*, and

_{M}*B*can be expressed as follows:where

*R*are calculated with the following form:where

*A*is the area of the pipe.

#### Boundary conditions

##### Pump

*H*and

_{PB}*H*correspond to the piezometric heads in the last and first sections of the suction and discharge pipes, respectively.

_{PA}*tdH*denotes the total dynamic head,

_{P}*Q*is the discharge of the pump,

_{P}*B*and

_{B}*B*are the pipeline constants

_{A}*B*of the suction and discharge pipes, respectively.

*tdH*and

_{P}*Q*can be obtained from the pump characteristic curve described by Suter (1966) by the linear interpolation. By substituting these values into Equation (12), the following equation is obtained:where

_{P}*α*=

*n*/

_{p}*n*and

_{r}*υ*=

*Q*/

_{p}*Q*.

_{r}*tdH*,

_{r}*n*, and

_{r}*Q*refer to the total dynamic head, rotational speed, and discharge under rated conditions, respectively.

_{r}*n*denotes the unknown rotational speed at the end of Δ

_{p}*t*.

*A*

_{0}and

*A*

_{1}denote the constants in the linear interpolation formula.

*β*

_{0}=

*T*

_{0}/

*T*.

_{r}*I*denotes the combined polar moment of inertia,

*α*

_{0}is the dimensionless speed at the beginning of Δ

*t*,

*T*

_{0}is the known torque at the beginning of Δ

*t*,

*T*is the rated torque on the pump, and

_{r}*T*is the unknown torque at the end of Δ

_{p}*t*.

*B*

_{0}and

*B*

_{1}are the constants contained in the linear interpolation formula.

##### Valve

*et al.*1993).where

*τ*is a dimensionless number that describes the discharge coefficient multiplied by the area of the valve opening. Δ

*H*

_{0}denotes the steady-state drop in the hydraulic grade line across the valve, which occurs at a flow rate of

*Q*

_{0}when

*τ*= 1. Δ

*H*denotes the steady-state drop with a flow rate of

*Q*, where

_{p}*Q*is the instantaneous discharge.

_{p}*B*

_{1}and

*B*

_{2}are the pipeline constants of pipes 1 and 2, respectively.

*C*

_{P}_{1}and

*C*

_{M}_{2}denote the known constants in the characteristic equations of pipes 1 and 2, respectively.

*C*is calculated as follows:

_{v}### Transient interference analysis of coupling relationships between multiple hydraulic facilities

*et al.*2021). One- and two-phase closure laws are the primary operating rules for valves in engineering applications (Zhao

*et al.*2019). Compared to the one-phase strategy, the two-phase strategy is advantageous in terms of limiting the flow and speed reversal of pumps, thereby reducing the amplitude of pressure pulsation and improving safety in the piping system. Additionally, in a series of hydraulic facilities and pipe systems, appropriately operating valves can effectively prevent excessive pressure and water levels during the successive shutdown process of pumps (Wan & Li 2016). Hence, the inputs

*P*comprise the following parameters: the time interval between the

*i*th and (

*i*

*+*1)th pumps shutdown (Δ

*T*

_{pi}_{−(i+1)}), total closing time (

*T*

_{vbp}_{2i}), inflection time (

*T*

_{vbp}_{1i}), and inflection opening (

*y*) of the valve closing law in the

_{vbpi}*i*th (

*i*= 1,2, … ,

*m*) pump-valve unit, the corresponding opening (

*y*) and required time to reach

_{rv_ji}*y*(

_{rv_ji}*T*) of the

_{rv_ji}*j*th (

*j*= 1,2, … ,

*t*) valve in the main pipeline when the

*i*th pump starts closing, as shown in Figure 2.

An ideal successive pump shutdown process minimizes the maximum water pressure (*HP*_{max}), highest water levels of the control structures (*Z _{e}*

_{max},

*e*

*=*1, 2,

*…*,

*s*) and total regulation time (

*T*

_{total}), as well as maximizing the minimum water pressure (

*HP*

_{min}) and lowest water levels (

*Z*

_{e}_{min}). Hence,

*HP*

_{max},

*HP*

_{min},

*Z*

_{e}_{max},

*Z*

_{e}_{min}, and

*T*

_{total}are selected as outputs

*O*, thereby ensuring the safety of the hydraulic transient process and quick feasibility of regulation.

### Multi-objective optimization

Aiming at the combined optimal operation of multiple hydraulic facilities based on NSGA-II and PNSGA-II, NSGA-II was utilized in a serial and parallel manner to drive the mathematical simulation model of the complex pressurized WDPs and calculate the corresponding real response functions (*N*_{pop}*×**O*) of the population (*N*_{pop}*×**P*). Then response functions were used to evaluate the implicit constraints and calculate objective functions.

However, when adopting the PMO-ASMO in the combined optimal operation model, a surrogate model focusing solely on the mathematical relationship between inputs *P* and outputs *O* is established to replace the simulation model. PNSGA-II was employed to reduce the computational burden of the optimization effectively. Further, an adaptive strategy was adopted to address the challenges of local optima and overfitting (Razavi *et al.* 2012), thereby enhancing the accuracy of optimization based on surrogate modeling. As shown in Figure 1, the mathematical simulation model is used to calculate the corresponding real response functions (*Z*_{0}*×**O*) of the initial samples (*Z*_{0}*×**P*), while the surrogate model is constructed using *Z*_{0}*×**O* and *Z*_{0}*×**P*. Then PNSGA-II runs on the surrogate model created in previous step to obtain the Pareto optimal points, matrix (*N*_{pop}*×**P*), where *N*_{pop} is the size of the population. A portion of the optimal points (20%) are chosen to solve the mathematical simulation model of complex pressurized WDPs. Further, the results were added to the initial data as (*Z*_{1}*×**P*) and (*Z*_{1}*×**O*) to update the surrogate model, wherein *Z*_{1} = *Z*_{0} + 20% *×**N*_{pop}. The surrogate model is iterated several times until the termination conditions are reached, which better describes the actual process of the physical model, thereby boosting the accuracy of the surrogate model and improving the global optimum.

#### Initial sampling

Latin hypercube (LH), Sobol’ low discrepancy sequence, Good Lattice Points (GLP), and crude Monte Carlo (MC) methods are the most commonly used sampling methods. The GLP with ranked Gram-Schmidt (RGS) decorrelation has been proven to perform optimally, since the computed discrepancy of GLP is small and RGS decorrelation can significantly improve the uniformity metrics of lattice designs (Gong *et al.* 2016b). Hence, GLP with RGS decorrelation was employed for initial sampling in this study, as described in more detail by Gong *et al.* (2016b). Furthermore, an appropriate initial sample size (*Z*_{0}) was set based on the specific situation.

*q*) and transmits them to the corresponding child thread of each computer core. A mathematical simulation model was used to calculate the corresponding response functions for the subsample set in each allocated child thread. Once the calculation of all the subsamples is complete, the corresponding response function matrix (

*Z*

_{0}

*×*

*O*) of the initial sample matrix (

*Z*

_{0}

*×*

*P*) is formed.

#### Surrogate model

The GPR, proposed by Rasmussen & Williams (2005), has been widely utilized to address nonlinear regression problems with small samples and high dimensionality. Hence, multiple surrogate models based on GPR were constructed for each response function in this study.

Moreover, the SCE-UA optimization method was chosen to determine the proper hyperparameter values in GPR (Gong *et al.* 2016a; Zhang *et al.* 2017).

#### Combined optimal operation model of multiple hydraulic facilities

Two objective functions, the minimization of the risk index for hydraulic transient processes and minimization of time consumption about regulation, were introduced in complex pressurized WDPs. This helps improve the safety and benefit of the combined operation of numerous hydraulic facilities during the successive shutdown process of the pumps. The aforementioned functions can be calculated based on the output *O*, whereas the decision variables are represented by the input *P*.

##### Objective functions

###### Minimization of risk index of hydraulic transient process **(***R*_{total})

*R*

_{total})

*et al.*2019), as follows:where

*HP*

_{s_}_{min}and

*HP*

_{s_}_{max}(

*HP*

_{s_}_{max}

*=*

*nHP*

_{w_}_{max}) represent the minimum and maximum allowable water pressure during the transient process, respectively, and

*HP*

_{w_}_{max}is the maximum steady-state pressure.

*n*is commonly set to 1.3, while

*HP*

_{s_}_{min}must not be less than −7.5 m to avoid cavitation.

*Z*,

_{es}*Z*

_{es}_{min}, and

*Z*

_{es}_{max}are the design, minimum, and maximum permissible water level of the

*e*th control structure, respectively.

*k*and

_{h}*k*represent the penalty factor for violating constraints of allowable water pressure and water levels of the

_{e}*e*th control structure, respectively.

###### Minimizing time consumption of regulation **(***T*_{total})

*T*

_{total})

*m*is the total number of pumps, Δ

*T*

_{pk−(k+}

_{1)}represents the time interval between the

*k*th and (

*k*

*+*1)th pumps shutdown,

*T*

_{vbp}

_{2m}is the total closing time for the valve closing law in the

*m*th pump-valve unit,

*T*

_{rv_jm}is the time required for the

*j*th valve on the main pipeline to reach corresponding opening (

*y*

_{rv_ji}) when the

*m*th pump starts to shut down.

##### Constraints

*O*derived from the mathematical simulation model or adaptive surrogate model.where

*t*

_{c_}_{min}and

*t*

_{c_}_{max}are the minimum and maximum allowable regulation times of the valves in the pump-valve units, respectively.

#### PNSGA-II

*et al.*2020), in which the fast-non-dominated sorting approach is used to greatly reduce computational complexity, by adopting the crowding distance to maintain the diversity of the population, and the elite strategy is introduced to make the sampling space enlarged and prevent the loss of the best individuals (Deb

*et al.*2002; Padhi & Mallick 2014; Hadibafekr

*et al.*2023). It has been widely applied in WDPs optimization operations. However, there is generally a high computational burden, even if a population with a large size is not difficult to implement. Moreover, the population of the NSGA-II utilizes the serial computing technique to assess the individual ability, thereby limiting the utilization efficiency of the computing resources owing to the ineffective utilization of parallel resources available in multicore machines. Consequently, when dealing with the combined operational problems of multiple hydraulic facilities in complex pressurized WDPs, the NSGA-II may require a substantial amount of time to complete the calculation (Peng

*et al.*2017). To avoid these defects, parallel computing was integrated into NSGA-II and the PNSGA-II was proposed (Liu

*et al.*2022). The Fork-Join framework, which is renowned for parallel computing (Niu

*et al.*2018), was used for parallel execution. The PNSGA-II workflow is illustrated in Figure 4, which mainly involves the following steps:

Step 1, a parent population, size *N*_{pop}, is randomly generated within the constraint ranges. Then it is arranged into *q* child threads, and the objective functions are calculated independently on every thread until completion. Furthermore, the computed results are returned to the master thread for being sorted into non-dominated ranks, while the crowding distances of each member are calculated.

Step 2, the parent population goes through the tournament selection, crossover, and mutation in NSGA-II, thereby resulting in a child population of size *N*_{pop}.

Step 3, *q* child threads are executed in parallel to obtain the objective functions of the child population. Further, all the results are returned to the master thread, and the parent and child population is merged. The 2*N*_{pop} solutions are sorted into non-dominated ranks. Similarly, the crowding distances are evaluated.

Step 4, *N*_{pop} members can be chosen as the parent population of next computational iteration from the 2*N*_{pop} population. Continue the aforementioned steps until the maximum iterations is reached.

## CASE STUDY

### Case description

*RV*

_{1}and

*RV*

_{2}are of a similar type.

*RV*

_{1}(or

*RV*

_{2}) and

*RV*

_{3}were closed by one-phase valve closure laws, where the time required to close their openings from 100 to 0% are 27 and 32 min, respectively. When pump 2 began shutting down, the valve on the branch pipe was closed using a one-phase valve closure law for 60 s. The characteristic parameters are presented in Table 1. Simultaneously considering the overall requirements for safe and efficient operation during the continuous shutdown of the three commonly used pumps is challenging. Consequently, designing a rationally combined operational scheme for multiple hydraulic facilities is important for enhancing the quality of the dynamic response.

Parameters . | Meaning . | Unit . | Value . | Parameters . | Meaning . | Unit . | Value . |
---|---|---|---|---|---|---|---|

Zur | Upstream reservoir level | m | 30.68 | Zdr | Downstream reservoir level | m | 41.17 |

L1 | Length of pipe 1 | m | 4,700 | D1 | Diameter of pipe 1 | m | 2 |

n1 | Manning roughness coefficient of pipe 1 | / | 0.015 | L2 | Length of pipe 2 | m | 29,400 |

D2 | Diameter of pipe 2 | m | 2.2 | n2 | Manning roughness coefficient of pipe 2 | / | 0.013 |

L3 | Length of pipe 3 | m | 2,150 | D3 | Diameter of pipe 3 | m | 2.6 |

n3 | Manning roughness coefficient of pipe 3 | / | 0.012 | L4 | Length of pipe 4 | m | 9,260 |

D4 | Diameter of pipe 4 | m | 2.2 | n4 | Manning roughness coefficient of pipe 4 | / | 0.012 |

L5 | Length of pipe 5 | m | 5,310 | D5 | Diameter of pipe 5 | m | 2.2 |

n5 | Manning roughness coefficient of pipe 5 | / | 0.013 | L6 | Length of pipe 6 | m | 5,490 |

D6 | Diameter of pipe 6 | m | 2 | n6 | Manning roughness coefficient of pipe 6 | / | 0.012 |

L7 | Length of pipe 7 | m | 4,050 | D7 | Diameter of pipe 7 | m | 2.2 |

n7 | Manning roughness coefficient of pipe 7 | / | 0.012 | L8 | Length of pipe 8 | m | 2,140 |

D8 | Diameter of pipe 8 | m | 2 | n8 | Manning roughness coefficient of pipe 8 | / | 0.012 |

Q1 | Discharge before the branch pipe | m^{3}/s | 5.5 | Q2 | Discharge after the branch pipe | m^{3}/s | 4.8 |

Qr | Pump rated discharge | m^{3}/s | 1.86 | Hr | Pump rated water head | m | 65.69 |

Nr | Pump rated speed | r/min | 750 | LZ1·WZ1 | Long·wide of the high-level water pool | m·m | 20 × 10 |

Z1smax | Allowable highest water level for the high-level water pool | m | 93.5 | Z1s | Design water level for the high-level water pool | m | 87.53 |

Z1smin | Allowable lowest allowable water level for the high-level water pool | m | 86.8 | LZ2·WZ2 | Long·wide of the regulating pool | m·m | 19.8 × 19.8 |

Z2smax | Allowable highest water level for the regulating pool | m | 65 | Z2s | Design water level for the regulating pool | m | 61 |

Z2smin | Allowable lowest water level for the regulating pool | m | 59.2 | n _{P10} | Initial ratio of operating speed to rated speed for pump 1 | – | 1 |

n _{p20} (n_{p30}) | Initial ratio of operating speed to rated speed for pumps 2 and 3 | – | 0.99 | τ(_{RV1}τ) _{RV2} | Initial opening degree of RV_{1} and RV_{2} | % | 73.87 |

τ _{RV3} | Initial opening degree of RV_{3} | % | 67.36 | – | – | – | – |

Parameters . | Meaning . | Unit . | Value . | Parameters . | Meaning . | Unit . | Value . |
---|---|---|---|---|---|---|---|

Zur | Upstream reservoir level | m | 30.68 | Zdr | Downstream reservoir level | m | 41.17 |

L1 | Length of pipe 1 | m | 4,700 | D1 | Diameter of pipe 1 | m | 2 |

n1 | Manning roughness coefficient of pipe 1 | / | 0.015 | L2 | Length of pipe 2 | m | 29,400 |

D2 | Diameter of pipe 2 | m | 2.2 | n2 | Manning roughness coefficient of pipe 2 | / | 0.013 |

L3 | Length of pipe 3 | m | 2,150 | D3 | Diameter of pipe 3 | m | 2.6 |

n3 | Manning roughness coefficient of pipe 3 | / | 0.012 | L4 | Length of pipe 4 | m | 9,260 |

D4 | Diameter of pipe 4 | m | 2.2 | n4 | Manning roughness coefficient of pipe 4 | / | 0.012 |

L5 | Length of pipe 5 | m | 5,310 | D5 | Diameter of pipe 5 | m | 2.2 |

n5 | Manning roughness coefficient of pipe 5 | / | 0.013 | L6 | Length of pipe 6 | m | 5,490 |

D6 | Diameter of pipe 6 | m | 2 | n6 | Manning roughness coefficient of pipe 6 | / | 0.012 |

L7 | Length of pipe 7 | m | 4,050 | D7 | Diameter of pipe 7 | m | 2.2 |

n7 | Manning roughness coefficient of pipe 7 | / | 0.012 | L8 | Length of pipe 8 | m | 2,140 |

D8 | Diameter of pipe 8 | m | 2 | n8 | Manning roughness coefficient of pipe 8 | / | 0.012 |

Q1 | Discharge before the branch pipe | m^{3}/s | 5.5 | Q2 | Discharge after the branch pipe | m^{3}/s | 4.8 |

Qr | Pump rated discharge | m^{3}/s | 1.86 | Hr | Pump rated water head | m | 65.69 |

Nr | Pump rated speed | r/min | 750 | LZ1·WZ1 | Long·wide of the high-level water pool | m·m | 20 × 10 |

Z1smax | Allowable highest water level for the high-level water pool | m | 93.5 | Z1s | Design water level for the high-level water pool | m | 87.53 |

Z1smin | Allowable lowest allowable water level for the high-level water pool | m | 86.8 | LZ2·WZ2 | Long·wide of the regulating pool | m·m | 19.8 × 19.8 |

Z2smax | Allowable highest water level for the regulating pool | m | 65 | Z2s | Design water level for the regulating pool | m | 61 |

Z2smin | Allowable lowest water level for the regulating pool | m | 59.2 | n _{P10} | Initial ratio of operating speed to rated speed for pump 1 | – | 1 |

n _{p20} (n_{p30}) | Initial ratio of operating speed to rated speed for pumps 2 and 3 | – | 0.99 | τ(_{RV1}τ) _{RV2} | Initial opening degree of RV_{1} and RV_{2} | % | 73.87 |

τ _{RV3} | Initial opening degree of RV_{3} | % | 67.36 | – | – | – | – |

### Results and discussion

As three frequently used pump-valve units are of the same type, *RV*_{1}, *RV*_{2}, and *RV*_{3} are closed at fixed speed, and the same type of valves, *RV*_{1} and *RV*_{2}, are installed in the same position and follow the same regulation, the decision variables are [*T _{vbp}*

_{21},

*T*

_{vbp}_{11},

*y*

_{vbp}_{1},

*ΔT*

_{p}_{1−2},

*ΔT*

_{p}_{2−3},

*y*

_{rv_}_{11},

*y*

_{rv_}_{12},

*y*

_{rv_}_{31},

*y*

_{rv_}_{32}]. The variation ranges of the constraints and decision variables are shown in Table 2.

Parameters . | Unit . | Values . | Parameters . | Unit . | Values . |
---|---|---|---|---|---|

t_{c_}_{min} | s | 10 | t_{c_}_{max} | s | 100 |

T_{vbp}_{21} | s | [10, 100] | y_{rv_}_{11} | % | [0, 73.87] |

T_{vbp}_{11} | s | [10,20] | y_{rv_}_{12} | % | [0, 73.87] |

Y_{vbp}_{1} | % | [10, 90] | y_{rv_}_{31} | % | [0, 67.36] |

ΔT_{p}_{1−2} | s | [390, 900] | y_{rv_}_{32} | % | [0, 67.36] |

ΔT_{p}_{2−3} | s | [390, 900] | – | – | – |

Parameters . | Unit . | Values . | Parameters . | Unit . | Values . |
---|---|---|---|---|---|

t_{c_}_{min} | s | 10 | t_{c_}_{max} | s | 100 |

T_{vbp}_{21} | s | [10, 100] | y_{rv_}_{11} | % | [0, 73.87] |

T_{vbp}_{11} | s | [10,20] | y_{rv_}_{12} | % | [0, 73.87] |

Y_{vbp}_{1} | % | [10, 90] | y_{rv_}_{31} | % | [0, 67.36] |

ΔT_{p}_{1−2} | s | [390, 900] | y_{rv_}_{32} | % | [0, 67.36] |

ΔT_{p}_{2−3} | s | [390, 900] | – | – | – |

The experimental setup of the algorithms is presented in Table 3. Considering the computational cost, the total number of runs (*N _{msm}*) for mathematical simulation model were chosen to be 400, 1000, and 1800. The setup of the NSGA-II was similar to that of PNSGA-II. For the PMO-ASMO, the maximum iteration (

*iter*) and population size (

*pop*) of the embedded PNSGA-II were both 100. Additionally, the size of adaptive sampling for every iteration is 100 × 20% = 20, wherein the ‘20%’ represents the specified resampling percentage (

*pct*).

*init*is the number of initial sampling.

N
. _{msm} | NSGA-II . | PNSGA-II . | PMO-ASMO . |
---|---|---|---|

400 | 20 pop × 20 iter | 20 pop × 20 iter | 200 init + (100 pop × 20% pct) × 10 iter |

1000 | 25 pop × 40 iter | 25 pop × 40 iter | 400 init + (100 pop × 20% pct) × 30 iter |

1800 | 30 pop × 60 iter | 30 pop × 60 iter | 800 init + (100 pop × 20% pct) × 50 iter |

N
. _{msm} | NSGA-II . | PNSGA-II . | PMO-ASMO . |
---|---|---|---|

400 | 20 pop × 20 iter | 20 pop × 20 iter | 200 init + (100 pop × 20% pct) × 10 iter |

1000 | 25 pop × 40 iter | 25 pop × 40 iter | 400 init + (100 pop × 20% pct) × 30 iter |

1800 | 30 pop × 60 iter | 30 pop × 60 iter | 800 init + (100 pop × 20% pct) × 50 iter |

The effectiveness of an algorithm could be evaluated based on the convergence and diversity of the Pareto optimal points relative to the theoretical Pareto frontier (Gong *et al.* 2016a). The number of evaluations for the function required for obtaining smooth Pareto optimal solution sets is regarded as the efficiency of an algorithm (Zhang *et al.* 2017). However, the theoretical Pareto frontier is unknown for the combined optimal operation of multiple hydraulic facilities in a true complex pressurized WDP. And performing numerous runs for complex pressurized WDPs simulations is not feasible owing to high computational costs. Therefore, in this study, the effectiveness of the PMO-ASMO was verified by comparing the performance of the Pareto solutions and the diversity of the Pareto solution set of the PMO-ASMO with that of the PNSGA-II and NSGA-II. The time consumed using PMO-ASMO, PNSGA-II, and NSGA-II in the combined optimal operation model of multiple hydraulic facilities was considered as the evaluation indicator of the efficiency.

#### Effectiveness

The purpose of the PNSGA-II algorithm is to effectively reduce the computational time by implementing multicore parallel computation while maintaining the same optimization effectiveness as the NSGA-II algorithm. Hence, with 400, 1000, and 1800 model runs, the Pareto optimal points obtained by PMO-ASMO was only compared to that of PNSGA-II.

*N*. As shown in Figure 7(a)–7(e), the number of Pareto solutions and feasible solutions of PMO-ASMO are less than those of PNSGA-II in the same number of simulation model runs. This is owing to the inevitable loss of some details during the search process, although the GPR surrogate model mimicked the effects of the simulation model. However, after 400 and 1000 simulations, the Pareto solutions of PMO-ASMO dominated those of PNSGA-II. While the Pareto solutions of the former partially dominated those of the latter in 1800 simulations, as observed during comparison. The reason for these results is that with the growth of the mathematical simulation model runs, the number of Pareto solutions also increases, and the distribution of the Pareto frontier solutions becomes more extensive, thereby encompassing more optimized Pareto frontier characteristics. Hence, PMO-ASMO, with limited Pareto solutions, cannot completely hold sway over the solutions of PNSGA-II when the number of simulation model runs increases. However, according to the box plot of the optimal objective functions from the 1800 model runs (Figure 8), the PMO-ASMO had the lowest median for the risk index objective function. Even though it had the highest median for regulation time using PMO-ASMO, the third quartile was the smallest, and the interquartile range the lowest. Although the Pareto solutions of PMO-ASMO may not entirely supersede those of PNSGA-II, they exhibit a better safety performance and are more appropriate for the actual requirements of engineering operations.

_{msm}Furthermore, because the GPR surrogate model exhibits the inevitable loss of some details during the search process, the proportions of feasible and Pareto solutions of PNSGA-II and PMO-ASMO to *N _{msm}* were studied to effectively determine

*N*, and the results are presented in Table 4. The numbers of feasible solutions of PNSGA-II with 400, 1000, and 1800 model runs were 149, 313, and 567, respectively, thereby accounting for 37.25, 31.30, and 31.50%, respectively, with an average of 33.35%. The number of Pareto solutions were 19, 21, and 43, thereby accounting for 4.75, 2.10, and 2.39%, respectively, with an average of 3.08%. For the feasible solutions with 400, 1000, and 1800 simulations in PMO-ASMO, the results were 19, 43, and 86, respectively, thereby representing 4.75, 4.30, and 4.77% of the total solutions, respectively, with an average of 4.61%. Further, the number of Pareto solutions were two, two, and seven, thereby accounting for 0.5, 0.2, and 0.39%, respectively, with an average of 0.36%. Hence, to ensure that feasible solutions of the PMO-ASMO contain a certain number,

_{msm}*N*is recommended to be not less than 220 in the combined optimal operation of multiple hydraulic facilities in complex pressurized WPDs.

_{msm}Algorithms . | N
. _{msm} | Numbers of feasible solutions . | Proportion of feasible solutions to N (%)
. _{msm} | Average proportion of feasible solutions to N (%)
. _{msm} | Numbers of Pareto frontier solutions . | Proportion of Pareto solutions to N (%)
. _{msm} | Average proportion of Pareto solutions to N (%)
. _{msm} |
---|---|---|---|---|---|---|---|

PNSGA-II | 400 | 149 | 37.25 | 33.35 | 19 | 4.75 | 3.08 |

1000 | 313 | 31.30 | 21 | 2.10 | |||

1800 | 567 | 31.50 | 43 | 2.39 | |||

PMO-ASMO | 400 | 19 | 4.75 | 4.61 | 2 | 0.5 | 0.36 |

1000 | 43 | 4.30 | 2 | 0.2 | |||

1800 | 86 | 4.77 | 7 | 0.39 |

Algorithms . | N
. _{msm} | Numbers of feasible solutions . | Proportion of feasible solutions to N (%)
. _{msm} | Average proportion of feasible solutions to N (%)
. _{msm} | Numbers of Pareto frontier solutions . | Proportion of Pareto solutions to N (%)
. _{msm} | Average proportion of Pareto solutions to N (%)
. _{msm} |
---|---|---|---|---|---|---|---|

PNSGA-II | 400 | 149 | 37.25 | 33.35 | 19 | 4.75 | 3.08 |

1000 | 313 | 31.30 | 21 | 2.10 | |||

1800 | 567 | 31.50 | 43 | 2.39 | |||

PMO-ASMO | 400 | 19 | 4.75 | 4.61 | 2 | 0.5 | 0.36 |

1000 | 43 | 4.30 | 2 | 0.2 | |||

1800 | 86 | 4.77 | 7 | 0.39 |

#### Efficiency

NSGA-II being time consuming in the combined optimal operation model of multiple hydraulic facilities, the efficiency of PMO-ASMO was compared to that of NSGA-II and PNSGA-II only for the scenario with 400 model runs. This comparison was made solely using PNSGA-II in 1000 and 1800 simulations.

In this study, it took approximately 980 s to execute one run of the mathematical simulation model for the case on a Core i9 (10 cores) 3.70 GHz processor (Intel, China). The time consumed using PNSGA-II and PMO-ASMO in the combined optimal operation model of multiple hydraulic facilities with 400, 1000, and 1800 model runs and the results after utilizing NSGA-II in the combined optimal operation model of multiple hydraulic facilities with 400 model runs are presented in Table 5. In the 400 model runs, PMO-ASMO exhibited a shorter execution time than the other two algorithms, whereas NSGA-II required the longest time. The time consumption of PNSGA-II with 10 cores in parallel was only 10.288% that of NSGA-II, thereby saving 103.625 h. This confirms the high efficiency of the PNSGA-II algorithm based on multicore parallelism. Compared to PNSGA-II, PMO-ASMO saved approximately 0.368 and 3.059 h in 400 and 1000 simulations, respectively. Nevertheless, PMO-ASMO did not demonstrate an advantage in terms of efficiency in the 1800 simulations as the computational time for PMO-ASMO was 62.224 h. However, it was only 54.577 h for PNSGA-II. Even though, PMO-ASMO can accelerate the searching process, which consistently acquires the best and stable Pareto frontier with fewer iterations in the combined operation model of multiple hydraulic facilities. Additionally, the computational time required by PMO-ASMO is several hours instead of days or even weeks as required by NSGA-II, thereby saving computing resources. Hence, PMO-ASMO can strengthen the applicability of the optimal operation of multiple hydraulic facilities in complex pressurized WDPs, and with continuous improvements in computing performances; it is also expected to be applied to the real-time optimal operation of multiple hydraulic facilities in complex pressurized WDPs.

N
. _{msm} | NSGA-II . | PNSGA-II . | PMO-ASMO . |
---|---|---|---|

400 | 415,832 (≈115.509 h) | 4,2783 (≈11.884 h) | 4,1456 (≈11.516 h) |

1000 | – | 121,913 (≈33.865 h) | 11,0902 (≈30.806 h) |

1800 | – | 196,476 (≈54.577 h) | 224,006 (≈62.224 h) |

N
. _{msm} | NSGA-II . | PNSGA-II . | PMO-ASMO . |
---|---|---|---|

400 | 415,832 (≈115.509 h) | 4,2783 (≈11.884 h) | 4,1456 (≈11.516 h) |

1000 | – | 121,913 (≈33.865 h) | 11,0902 (≈30.806 h) |

1800 | – | 196,476 (≈54.577 h) | 224,006 (≈62.224 h) |

#### Analysis of solutions

Considering the efficiency requirements in an actual operation, the solutions corresponding to the optimal value of each objective function in the Pareto fronts obtained PMO-ASMO and PNSGA-II with 400 simulations (Figure 9) were selected as the optimal operation schemes. The schemes corresponding to the minimum values of the response and objective functions are presented in Table 6.

. | PNSGA-II (N = 400)_{msm}. | PMO-ASMO (N = 400)_{msm}. | |||
---|---|---|---|---|---|

R: min and _{total}T: max
. _{total} | R: max and _{total}T: min
. _{total} | R: min and _{total}T: max
. _{total} | R: max and _{total}T: min
. _{total} | ||

Decision variables | T_{vbp}_{11} (s) | 11 | 11 | 16 | 18 |

y_{vbp}_{1} | 18.77% | 18.77% | 19.87% | 19.63% | |

T_{vbp}_{21} (s) | 79 | 79 | 74 | 73 | |

ΔT_{p}_{1−2} (s) | 730 | 680 | 570 | 510 | |

ΔT_{p}_{2−3} (s) | 640 | 620 | 410 | 600 | |

y_{rv}_{_11} | 58.32% | 54.52% | 58.24% | 57.08% | |

y_{rv}_{_12} | 18.81% | 16.24% | 33.48% | 20.04% | |

y_{rv}_{_31} | 29.77% | 31.94% | 38.14% | 40.80% | |

y_{rv}_{_32} | 12.57% | 10.89% | 23.61% | 9.58% | |

Response functions | HP_{max} (m) | 86.98 | 88.24 | 86.44 | 87.07 |

HP_{min} (m) | 1.20 | 1.20 | 1.27 | 1.29 | |

Z_{1max} (m) | 90.06 | 91.24 | 88.08 | 90.11 | |

Z_{1min} (m) | 86.81 | 86.80 | 86.87 | 86.89 | |

Z_{2max} (m) | 60.99 | 60.91 | 60.91 | 60.91 | |

Z_{2min} (m) | 59.67 | 59.27 | 59.46 | 59.24 | |

T_{total} (s) | 1,675 | 1,563 | 1,522 | 1,435 | |

Objective functions | R_{total} | 1.37 | 1.64 | 1.11 | 1.45 |

T_{total} (s) | 1,675 | 1,563 | 1,522 | 1,435 |

. | PNSGA-II (N = 400)_{msm}. | PMO-ASMO (N = 400)_{msm}. | |||
---|---|---|---|---|---|

R: min and _{total}T: max
. _{total} | R: max and _{total}T: min
. _{total} | R: min and _{total}T: max
. _{total} | R: max and _{total}T: min
. _{total} | ||

Decision variables | T_{vbp}_{11} (s) | 11 | 11 | 16 | 18 |

y_{vbp}_{1} | 18.77% | 18.77% | 19.87% | 19.63% | |

T_{vbp}_{21} (s) | 79 | 79 | 74 | 73 | |

ΔT_{p}_{1−2} (s) | 730 | 680 | 570 | 510 | |

ΔT_{p}_{2−3} (s) | 640 | 620 | 410 | 600 | |

y_{rv}_{_11} | 58.32% | 54.52% | 58.24% | 57.08% | |

y_{rv}_{_12} | 18.81% | 16.24% | 33.48% | 20.04% | |

y_{rv}_{_31} | 29.77% | 31.94% | 38.14% | 40.80% | |

y_{rv}_{_32} | 12.57% | 10.89% | 23.61% | 9.58% | |

Response functions | HP_{max} (m) | 86.98 | 88.24 | 86.44 | 87.07 |

HP_{min} (m) | 1.20 | 1.20 | 1.27 | 1.29 | |

Z_{1max} (m) | 90.06 | 91.24 | 88.08 | 90.11 | |

Z_{1min} (m) | 86.81 | 86.80 | 86.87 | 86.89 | |

Z_{2max} (m) | 60.99 | 60.91 | 60.91 | 60.91 | |

Z_{2min} (m) | 59.67 | 59.27 | 59.46 | 59.24 | |

T_{total} (s) | 1,675 | 1,563 | 1,522 | 1,435 | |

Objective functions | R_{total} | 1.37 | 1.64 | 1.11 | 1.45 |

T_{total} (s) | 1,675 | 1,563 | 1,522 | 1,435 |

As shown in Table 6, the lowest risk index scheme computed by PMO-ASMO had a risk index of 1.11 and a duration of the regulation of 1,522 s, which were reduced by 18.98 and 9.13%, respectively, compared to those obtained by PNSGA-II. For the shortest regulation time scheme calculated using PMO-ASMO, the risk index was 1.45 and duration of regulation 1,435 s, which were reduced by 11.59 and 8.19%, respectively, compared to the scheme calculated using PNSGA-II. The reason is that the response function values in the two schemes calculated using PMO-ASMO are improved except for only a slight decrease in the lowest water level of the regulating pool. In particular, the water level fluctuations relative to the designed water level decreased by 60.92 and 27.48% in the high-level water pool for the lowest risk index scheme and the shortest regulation time scheme, respectively.

The total shut time of the valve closing law in the pump-valve unit is shorter relative to the total regulation time, as well as the fact that the corresponding openings of *RV*_{3} when each pump starts closing has less impact on the high-level water pool due to the presence of regulating pool. The main factors affecting the change of water level in the high-level water pool are the time intervals between pumps shutdown, and the corresponding openings of *RV*_{1} (or *RV*_{2}). Therefore, the fundamental reason why the schemes calculated by PMO-ASMO are all better than obtained by PNSGA-II, is that the difference in water discharge into and out of the high-level water pool is minimized by reasonably adjusting the corresponding openings of *RV*_{1} (or *RV*_{2}) under shorter pumps shutdown intervals, which in turn reduces the water level fluctuations. Moreover, the relationship among the closing law of the valve in the pump-valve unit, the corresponding openings of *RV*_{1} (or *RV*_{2}) and *RV*_{3}, and the time intervals between pumps shutdown is also coordinated at the same time, so that the system pressure extremes and the water level fluctuations of regulating pool do not get worse, or even improve.

## CONCLUSIONS

This study focused on achieving a combined optimal operation of multiple hydraulic facilities for complex pressurized WDPs. The combined optimal operation model of the multiple hydraulic facilities coupled with mathematical simulation models of complex pressurized WDPs is established based on the successive shutdown process of pumps, and the PMO-ASMO algorithm that combines adaptive surrogate models and parallel computing is provided. Using the complex pressurized diversion section of the Jiaodong Water Diversion Project as an example, the results of the PMO-ASMO, NSGA-II, and PNSGA-II algorithms were compared. The findings are as follows:

(1) The optimization effectiveness of PMO-ASMO outperformed PNSGA-II not only by using similar number of mathematical simulations, but also by using fewer model runs.

(2) The time consumption of PMO-ASMO was only 9.97% of that acquired by NSGA-II, which was comparable to that of PNSGA-II, in the case of 10-core parallelism.

(3) The PMO-ASMO outperforms NSGA-II and PNSGA-II with the same number of mathematical simulations owing to the following reasons: (a) the adaptive surrogate model is employed to simulate the complex pressurized WDPs to alleviate the computational burden, while (b) the PNSGA-II algorithm can realize multicore parallel computation.

Therefore, a high-efficient and economical framework is provided for obtaining combined optimal operation schemes of multiple hydraulic facilities with the different operational conditions in WDPs. However, there exist some limitations caused by insufficient diversity in the algorithm settings. Further research is required to explore the combined operation of multiple hydraulic facilities under other operating conditions.

## ACKNOWLEDGEMENTS

This work was partly supported by the Basic Research Programs of Shanxi Province (grant nos. 20210302124645, 202203021222112), the National Key R&D Program of China (grant no. 2021YFC3001000), and the Open Research Fund of Henan Key Laboratory of Water Resources Conservation and Intensive Utilization in the Yellow River Basin (grant no. HAKF202104).

## DATA AVAILABILITY STATEMENT

All relevant data are included in the paper or its Supplementary Information.

## CONFLICT OF INTEREST

The authors declare there is no conflict.