Short-Term Multi-Objective Optimal Operation of Reservoirs to Maximize the Benefits of Hydropower and Navigation

Abstract

Traditional reservoir operation mainly focuses on economic benefits, while ignoring the impacts on navigation. Thus, the economic operation of reservoirs considering navigational demands is of great significance for improving benefits. A navigation capacity evaluation method (NCEM), which evaluates the navigation capacity considering the influence of flow velocity and water level variation on navigation, is proposed to more effectively evaluate the navigation capacity. Based on two-dimensional hydrodynamic numerical simulation, the NCEM accurately calculates the navigation capacity according to detailed flow velocity and water level changes. In addition, a short-term multi-objective optimal operation model considering the upstream and downstream navigation and power generation is established. Then, the Strength Pareto Evolutionary Algorithm (SPEA2) is used to solve the model. To verify the rationality of the method and model, they are applied to the case study of the Xiangjiaba reservoir. The results demonstrate that the method and model can not only provide a series of operation schemes for decision makers of reservoirs, but also direct the ship to pass safely through the approach channel, implying a certain practical value and significance as a reference for the short-term optimal operation of reservoirs in the future.

1. Introduction

Hydropower is an important clean and renewable energy source, and increasing numbers of hydropower stations have been built in recent years to meet economic needs worldwide [1,2,3,4,5,6]. In order to explore resources more efficiently, many scholars have carried out studies on the optimal operation of reservoirs [7,8,9]. Furthermore, some useful optimal models and algorithms which comprehensively consider the benefits of reservoirs are proposed [10,11,12,13]. However, there are relatively few studies on the operation of reservoir considering navigation demands.

During the daily operation of the reservoir, scheduling schemes are developed to take into account power generation needs [14], peak load regulation [15], and storage water demands etc., which will result in daily cyclical changes in the discharge volume and water level [16]. The fluctuation is transmitted downstream, causing a complex unsteady flow [17]. Moreover, due to the special topography of some rivers, the large discharge volume of the reservoir can cause a large surface velocity and reflux in the local area of the river reaches [18]. These phenomena seriously hinder the navigation of ships [19]. Consequently, in order to ensure the safe and stable operation of the reservoir, it is necessary to implement an evaluation of the navigation capacity of the channel, laying a foundation for a comprehensive optimization of the operation of the reservoir [20].

Defining and quantifying navigation capacity is the premise for the short-term optimal operation of reservoirs to meet the requirements of navigation, and it is also the basis for assessing the impact of operation of the reservoir on navigation. There are many studies on the operation mode of ship locks in order to increase navigation benefits [21,22,23,24], defining and quantifying the navigation capacity to a certain extent so that the navigation objective can be integrated into the optimization model. Ji et al. [25,26,27] used the weighted total tardiness of the ships as the navigation objective to solve the problem of ship lock scheduling, increasing the interests of shippers and owners of the lock. To maximize navigation benefits, the Non-Dominated Sorting Genetic Algorithm III (NSGA-III) was used to optimize the operation mode of the ship lock and the order of the ship passing through the ship lock. However, this method mainly focuses on the navigation efficiency and the operation mode of the ship lock. It does not highlight the operation effect of the reservoir and lacks a comprehensive consideration of the multiple benefits of the reservoir.

Taking the hydraulic factors of the reservoir and river reaches as indicators of navigation capacity, the navigation objective and other objectives (power generation, flood control, ecological environment, etc.) of the reservoir can be combined to form a multi-objective optimization model. The influencing factors of the navigation are complex, and parameters such as flow velocity, water surface gradient, variation of daily water level, and variation of hourly water level are usually used as the evaluation indicators at present. Ackermann et al. [28] used the sum of water level variability and flow variability with different weighting factors as the navigation objective based on the result of the one-dimensional hydrodynamic simulation. Wang and Zhang [29] proposed a single objective optimization model considering navigation demands (using the maximum amount of power generation as the objective), and took flow ramping as the constraint. Ma et al. [30] used the comprehensive influence of multiple hydraulic factors for the evaluation of navigability, defining the navigation capacity as zero (not satisfying the navigation requirements) or one (satisfying the navigation requirements) when the navigation capacity could not be evaluated in detail. Liu et al. [31] used the guarantee rate of navigation to quantify the navigation objective, which is related to the discharge volume. This method combines well with the flood control objective and navigation objective, meets the demands for the multiple benefits of reservoir, and achieves good results.

In summary, the above studies analyze the influence of reservoir scheduling on navigation, and use intelligent optimization algorithms to solve optimization problems. In addition, optimal operation schemes are achieved for both single-objective and multi-objective models. Nevertheless, the above methods and models at least have one of the following shortcomings. Firstly, the navigation objective may not be defined and quantified, and the navigation conditions are treated as constraints rather than objectives in the model. The obtained optimization results thus do not reflect the navigation benefits. Secondly, the quantization of the navigation objective may be too rough to consider the impact of some hydraulic factors on navigation accurately. Thirdly, the evaluation of the navigation capacity may not consider the impact of flow velocity on navigation, and may not take into account the combined effects of various hydraulic factors (flow rate, water level, etc.) on navigation. Fourthly, a single-objective model may be used to analyze the comprehensive benefits of navigation and power generation, which has great limitations in exploring the relationship between multiple objectives.

To this end, this paper proposes a navigation capacity evaluation method which evaluates the navigation capacity considering the comprehensive influence of flow velocity and water level variation on navigation. Moreover, a short-term multi-objective optimal operation model considering the upstream and downstream navigation and power generation is established.

The rest of this paper is organized as follows. Section 2 firstly analyzes the impact of the flow velocity and water level on navigation, and a navigation capacity evaluation method (NCEM) is proposed for defining and quantifying the navigation objective. In Section 3, the multi-objective operation model of reservoir is established. Section 4 briefly shows the algorithm and flowchart of the multi-objective model adopted in this study. Section 5 presents the case study of the Xiangjiaba reservoir, with calculation results and discussions. Finally, the conclusions are given in Section 6.

2. Navigation Capacity Evaluation

The navigation capacity evaluation method (NCEM) which evaluates the navigation capacity considering the comprehensive influence of flow velocity and water level variation on navigation of ships, is proposed in this paper. The navigation capability in NCEM is composed of the navigation capacity considering the flow velocity ($NCv$) and the navigation capacity considering the water level variation ($NCl$) according to a certain weighting factor. Based on the hydrodynamic simulation, the flow field and water level variation in the upstream and downstream river reaches can be calculated and accurately estimated.

2.1. Hydrodynamic Model

Based on the constant fluid density, incompressible characteristics, and Navier–Stokes formulas, a two-dimensional hydrodynamic model is used to analyze the unsteady flow [32]. The corresponding continuity formula is as follows:

$$\frac{\partial h}{\partial t}+\frac{\partial hu}{\partial x}+\frac{\partial hv}{\partial y}=S$$

(1)

The momentum formulas in the $x$ and $y$ directions are as follows:

$$\frac{\partial u}{\partial t}+u\frac{\partial u}{\partial x}+v\frac{\partial u}{\partial y}=-g\frac{\partial \eta }{\partial x}-g{S}_x+\frac{1}{h}\frac{\partial }{\partial x}\left(v_{t}h\frac{\partial u}{\partial x}\right)+\frac{1}{h}\frac{\partial }{\partial y}\left(v_{t}h\frac{\partial u}{\partial y}\right)$$

(2)

$$\frac{\partial v}{\partial t}+u\frac{\partial v}{\partial x}+v\frac{\partial v}{\partial y}=-g\frac{\partial \eta }{\partial y}-g{S}_y+\frac{1}{h}\frac{\partial }{\partial x}\left(v_{t}h\frac{\partial v}{\partial x}\right)+\frac{1}{h}\frac{\partial }{\partial y}\left(v_{t}h\frac{\partial v}{\partial y}\right)$$

(3)

where $u\left(x,y,t\right)$ and $v\left(x,y,t\right)$ are the velocity components in the horizontal $x$ and $y$ directions, respectively; $S$ denotes source terms involving bed slope and friction; $h$ denotes the total depth of water measured from bed; $t$ denotes the time; $\eta (x,y,t)$ represents the water surface elevation; $g$ is the gravitational acceleration; and $v_t$ represents the coefficients of horizontal eddy viscosity in the momentum formula. $S_x$ and $S_y$ are bed resistance in the $x$ and $y$ directions, respectively, which can be expressed as follows:

$$S_x=\frac{n^{2}u\sqrt{u^2+v^2}}{h^{4/3}}$$

(4)

$$S_y=\frac{n^{2}v\sqrt{u^2+v^2}}{h^{4/3}}$$

(5)

where $n$ is the Manning’s roughness coefficient.

In this paper, the software Mike21 FM is used for modeling and numerical simulation. In order to make the mathematical model adapt to the complex boundary of natural rivers with high computational efficiency, a two-dimensional shallow water formula based on the unstructured grid is carried out in this study. The continuity formula is discretized by the finite volume method [33]. The two-stage TVD Runge–Kutta method is developed to solve the shallow water formulas in the model [34,35].

2.2. Analysis of the Flow Velocity

During the operation of the reservoir, excessive transverse velocity will be generated in the entrance area of approach channel, thus hindering navigation. Therefore, it is of importance to determine the dangerous navigation area in the entrance area of the approach channel under different operational conditions of the reservoir. The determination of the dangerous area will lay a solid foundation for the evaluation of navigation capacity considering the flow velocity ($NCv$).

2.2.1. The Dangerous Navigation Area of Channel

In the actual navigation process, as long as the width of the area where the transverse velocity exceeds the limit value is no more than two-thirds of the length of the ship ($d_{ship}$), the ship can pass through this area safely.

In the transverse velocity distribution map obtained from hydrodynamic simulation, $Dc$ is defined as the width of the approach channel, and the dangerous width of channel ($Dd$) is expressed as the maximum width of the dangerous navigation area. Moreover, the dangerous navigation area is defined as the region where not only does the transverse velocity exceed the maximum limit velocity but also the length of the area is not less than two-thirds of $d_{ship}$ [36]. A graphical interpretation of the dangerous navigation area and $Dd$ are shown in Figure 1.

2.2.2. The Navigation Capacity Considering the Flow Velocity ($NCv$)

The value of $Dd$ is determined by the value of $d_{ship}$ and the operation condition of the reservoir. In addition, the larger the value of $Dd$, the worse the navigation environment of the channel, and the smaller the navigation capacity considering the flow velocity. Thus, the $NCv$ can be expressed by the following formulas:

$$NC{v}_t=1-\frac{1}{Dc}Dd(Q{x}_t,d_{ship})$$

(6)

$$NCv=\frac{1}{Tn}{\displaystyle \sum _{t=1}^{Tn}\left[1-\frac{1}{Dc}Dd(Q{x}_t,d_{ship})\right]}$$

(7)

where $NC{v}_t$ represents the navigation capacity considering the flow velocity in the $t$-th stage; $NCv$ represents the average value of $NC{v}_t$ in $Tn$ time periods; $Tn$ is the total number of periods in the navigation period; $Dc$ is the width of the approach channel; and $Q{x}_t$ is the discharge volume in the $t$-th stage.

2.3. Analysis of the Variation of Water Level

When ensuring that the influence of the flow velocity in the river channel or the approach channel is not significant, it should be confirmed that the water level is stable and the amplitude should not be too large for navigation. The hourly variation of the water level is an important measure of navigation capacity. Furthermore, the smaller the value of the water level variation, the greater the navigational capacity.

Through the solution of the hydrodynamic model, it can be clearly seen that the water level at different discrete cross-sections changes with the operation of the reservoir, and then the navigation capacity considering the water level variation ($NCl$) can be obtained. The expression of $NCl$ is as follows:

$$NC{l}_t=\frac{1}{CS}{\displaystyle \sum _{j=1}^{CS}(1-\left|\frac{Z_{t,j}-Z_{t-1,j}}{\Delta Z_{h,\lim \mathrm{t}}}\right|)}$$

(8)

$$NCl=\frac{1}{CS}\frac{1}{\left(Tn-1\right)}{\displaystyle \sum _{j=1}^{CS}{\displaystyle \sum _{t=2}^{Tn}(1-\left|\frac{Z_{t,j}-Z_{t-1,j}}{\Delta Z_{h,\lim \mathrm{t}}}\right|)}}$$

(9)

where $NC{l}_t$ represents the navigation capacity considering the water level variation in the $t$-th stage; $NCl$ represents the average value of $NC{l}_t$ in $Tn$ time periods; $CS$ is the number of the discrete cross-sections of the river reach; $Z_{t,j}$ is the water level of the $j$-th cross-section in the $t$-th stage; and $\Delta Z_{h,\lim \mathrm{t}}$ is the maximum limit variation per hour of water level.

3. Multi-Objective Reservoir Operation Model

The short-term multi-objective optimal operation model of reservoir considers economic and navigation benefits in this paper. The model consists of three objectives: the economic objective, the upstream navigation capacity objective, and the downstream navigation capacity objective. The decision variables in this paper are set to the upstream water level ($\left\{Z_1^{up},Z_2^{up},Z_3^{up},\dots ,Z_{T+1}^{up}\right\}$) of the reservoir.

The power generation of the reservoir is directly related to the flow used to generate electricity, and has an indirect relationship with the discharge volume. The navigation capacity is directly related to the discharge volume. Therefore, the discharge volume is used as the link to analyze the influence of short-term reservoir operation on power generation and navigation.

3.1. Objective Function

3.1.1. Economic Objective: Maximizing the Daily Total Power Generation

Maximizing the total daily power generation is a necessary objective to guarantee economic benefits [37]. For short-term optimal operation problem, the formula for economic objective can be written as follows:

$$\max E={\displaystyle \sum _{t=1}^T{N}_t\left(Q_t,H_t\right)}\times \Delta t,t\in \left[1,T\right]$$

(10)

where $T$ is the total number of periods in the operation period; $\Delta t$ is the length of time in each period; $E$ is the total power generation; $N_t$ is the output of hydroelectric generators in the $t$-th stage; and $Q_t$ and $H_t$ are the power discharge and net head in the $t$-th stage, respectively.

Moreover, $N_t$ depends on $Q_t$, $H_t$, and the energy loss, and the energy loss is determined by the efficiency of the hydroelectric generator. The formula can be expressed as follows:

$$N_t=K{Q}_t{H}_t$$

(11)

$$H_t=Z_t^{up}-Z_t^{down}-H_{lose}$$

(12)

where $Z_t^{up}$ and $Z_t^{down}$ are the water level in the upstream and downstream of the reservoir in the $t$-th stage, respectively; and $H_{lose}$ is the water head loss.

3.1.2. Navigation Objective: Maximizing the Navigation Capacity

To maximize the navigation capacity is an important objective, which is beneficial to navigation. Besides, the navigation capacity is a comprehensive evaluation indicator, taking into account the influence of flow velocity and water level variation on the ship according to Section 2. The formula to calculate the navigation objective is defined as follows:

$$\max NC={\omega }_v\times NCv+{\omega }_l\times NCl$$

(13)

$${\omega }_v+{\omega }_l=1$$

(14)

where $NC$ represents the navigation capacity; and ${\omega }_v$ and ${\omega }_l$ are weighting factors. According to Formulas (7) and (9), Formula (13) can also be written as the following form:

$$NC={\omega }_v\times \frac{1}{Tn}{\displaystyle \sum _{t=1}^{Tn}\left[1-\frac{1}{Dc}Dd(Q{x}_t,d_{ship})\right]}+{\omega }_l\times \frac{1}{CS}\frac{1}{\left(Tn-1\right)}{\displaystyle \sum _{j=1}^{CS}{\displaystyle \sum _{t=2}^{Tn}(1-\left|\frac{Z_{t,j}-Z_{t-1,j}}{\Delta Z_{h,\lim \mathrm{t}}}\right|)}},t\in \left[1,T_n\right]$$

(15)

3.2. Constraints

Due to the operation of reservoir and the demands of navigation, the optimization of three objectives is restricted due to the following constraints [38].

(1) Water volume balance

$$V_t=V_{t-1}+\left(I_t-Q{x}_t\right)\times \Delta t,t\in \left[1,T\right]$$

(16)

where $V_t$ is the storage volume of the reservoir in the $t$-th stage and $I_t$ is the inflow in the $t$-th stage.

(2) Constraints on the reservoir volume

$$V_t^{\min }\le V_t\le V_t^{\max }$$

(17)

where $V_t^{\min }$ is the minimum limit of the $V_t$ in the $t$-th stage and $V_t^{\max }$ is the maximum limit of the $V_t$ in the $t$-th stage.

(3) Water head constraint of hydropower station:

$$H_{\min }\le H_t\le H_{\max }$$

(18)

where $H_{\min }$ and $H_{\max }$ are the minimum and maximum limit of work head level in the $t$-th stage respectively, which depends on the hydroelectric generator.

(4) Power generation constraint:

$$N_{\min }\le N_t\le N_{\max }$$

(19)

where $N_{\min }$ and $N_{\max }$ are the lower and upper limit of $N_t$ in the $t$-th stage, respectively.

(5) Navigable discharge volume constraint

$$Q{n}_t^{\min }\le Q{x}_t\le Q{n}_t^{\max },t\in \left[1,T_n\right]$$

(20)

where $Q{n}_t^{\min }$ and $Q{n}_t^{\max }$ are minimum and maximum navigable discharge volume in the $t$-th stage, respectively.

(6) Boundary conditions constraint

$$\{\begin{array}{l}{Z}_1=Z^{\mathrm{beg}}\\ Z_{T+1}=Z^{\mathrm{end}}\end{array}$$

(21)

where $Z^{\mathrm{beg}}$ is the upstream water level in the beginning of the scheduling period; $Z^{\mathrm{end}}$ is the upstream water level in the end of the scheduling period.

4. Methodology

Studying the relationship between power generation, the upstream navigation capacity, and the downstream navigation capacity is a multi-objective optimization problem (MOP) [39,40]. There are several multi-objective evolutionary algorithms to solve the complicated multi-objective problem, such as the Strength Pareto Evolutionary Algorithm (SPEA2) [41]; the Non-Dominated Sorting Genetic Algorithm II (NSGA-II) [42]; the Multi-Objective Bee Colony Optimization algorithm (EMOBCO) [43]; and the Multi-Objective Artificial Sheep Algorithm (MOASA) [44], etc. Besides, maintaining a good distribution of solutions is a challenge for the optimization algorithms. Therefore, the SPEA2, which is a multi-objective evolutionary algorithm based on elitist strategy, is used to solve the MOP problem with three objectives in this study. The SPEA2 employs an enhanced fitness assignment strategy for archive truncation and density-based selection in order to obtain the Pareto frontier with uniform distribution. In the algorithm, the archive is used to store all non-dominated solutions generated in the optimization process, and the advantage of the archive is to obtain the Pareto front with better distribution. To start with an initial population and an empty archive, all non-dominated individuals are copied to the archive, and the dominated population members are removed from the archive under the update strategy of the archive set for each iteration. The update strategy of the archive set is as follows:

(1)

When there are no individuals in the external archive set, the non-dominated solutions generated in the iteration are stored directly in the archive set;

(2)

If the newly generated individuals are dominated by the individuals in the archive set, the new individuals will be deleted. Conversely, individuals can be deleted from the original archive set and new individuals added to the archive set;

(3)

When the number of individuals in the archive set reaches the preset maximum capacity, the individuals with smaller crowding distance are deleted.

Therefore, the Pareto solution set is stored in the archive, and the archive set is output as the optimization result finally.

Coupled with the above SPEA2 and NCEM, the flowchart of the multi-objective model is shown in Figure 2.

The calculation steps are as follows:

Step 1: Establish the hydrodynamic model

Two two-dimensional hydrodynamic models (upstream and downstream river reaches) are established according to Section 2.1. The flow profile in the unsteady flow condition is used as the upstream boundary, and the water level profile is used as the downstream boundary in both models [17]. Based on the results of the numerical simulation, the value of $Dd$ and the cross-sectional water level corresponding to different operation conditions are calculated. Therefore, the value of $NCv$ and $NCl$ under different operation conditions can be determined according to Formulas (6) and (8).

Step 2: Determine the weighting factors

${\omega }_v$ and ${\omega }_l$ are weighting factors that determine the effect of flow velocity and water level variation on navigation capacity in Formula (13).

The value of ${\omega }_v$ and ${\omega }_l$ can be determined according to Step 1. If there is very large flow velocity in the river reach during the operation of reservoir, then the flow velocity is considered to be the main obstacle to navigation, and the value of ${\omega }_v$ is set greater than ${\omega }_l$. If there is no large flow velocity, the water level variation is the main factor affecting navigation, and the value of ${\omega }_l$ is greater than ${\omega }_v$.

Step 3: Determine the navigation objective functions

According to Formula (13), the parameters (${\omega }_v$) and (${\omega }_l$) of the upstream and downstream navigation objective functions are determined in Step 2. In the calculation of each navigation objective function of the following algorithm, the value of $NCv$ and $NCl$ under different operation conditions calculated in Step 1 will be imposed in the navigation objective function to determine the specific value of navigation capacity.

Step 4: Initialization

A series of individuals which represent the upstream water level $\left\{Z_1^{up},Z_2^{up},Z_3^{up},\dots ,Z_{T+1}^{up}\right\}$ are generated, and these individuals are randomly generated according to maximum and minimum limits in the algorithm. Besides, the parameters of SPEA2 are initialized, including the population size, the maximum number of generations, the size of the archive, the crossover rate, the mutation rate, etc.

Step 5: Solution algorithm

In this step, the multi-objective problem is solved by SPEA2. Among them, the expression form of the navigation objective function is determined by the result calculated in Step 3. The value of $NCv$ and $NCl$ in different operation conditions can be obtained by the numerical simulation in Step 1, providing a basis for the calculation of each navigation objective function.

Step 6: The output

The archive set is output as the optimization result if the termination condition is met, and the Pareto solution set is obtained finally.

5. Case Study

5.1. Case Description

The Xiangjiaba (XJB) reservoir is located at the juncture of Shuifu city and Yibin city in China, which is the last level reservoir of Jinsha River hydropower base, as shown in Figure 3a,b. With a multi-year average runoff of 3810 m³/s, the XJB reservoir can control a basin area of 458,800 square kilometers, accounting for 97% of the Jinsha River basin area. In addition, the total storage capacity of the XJB reservoir is 5.163 billion cubic meters, and the regulating storage volume is 900 million cubic meters. The main parameters of XJB reservoir is shown in

Table 1. Basic parameters of the XJB reservoir.

Characteristics of Water Level (m)		Installed Capacity/MW	Minimum Power/MW	$\mathit{Q}{\mathit{n}}_{\mathbf{max}}({\mathbf{m}}^3/\mathbf{s})$	$\mathit{Q}{\mathit{n}}_{\mathbf{min}}({\mathbf{m}}^3/\mathbf{s})$	Flow Velocity (m/s)		$\mathbf{\Delta }{\mathit{Z}}_{\mathit{h},\mathbf{lim}\mathit{t}}(\mathbf{m}/\mathbf{h})$
Dead Water Level	The Normal Height					umax1 (m/s)	vmax2 (m/s)
370	380	6400	1800	12,000	1200	2.0	0.3	1.5

¹u_max is the maximum limit longitudinal velocity; ² v_max is the maximum limit transverse velocity.

.

The XJB reservoir not only has the function of generating electricity, but also regulating navigation. It is the third largest hydropower station in China and the fifth largest in the world, for which the total installed capacity for power generation is 6 × 10⁶ kW, and the annual average power generation is 3.013 × 10¹⁰ kWh. Furthermore, it is the only dam in the Jinsha River hydropower base to build a ship lift with a very efficient ship navigation rate. Its scale is equivalent to that of the Three Gorges reservoir, and it takes only 15 min for a ship of one kiloton to cross the dam. Compared with the average dam time of the Three Gorges reservoir of 5 h, the navigation efficiency of XJB reservoir is quite high.

However, due to the operation of reservoir, the navigation of the upstream and downstream reaches of the reservoir is seriously affected by the unsteady flow and the navigational is impeded. Besides, the difficulty of upstream and downstream navigation of the reservoir is both caused by the water flow at the entrance area of the approach channel. When the ship passes through the entrance area, care should be taken to avoid the influence of the large transverse velocity near the dike on the navigation. Therefore, it is of great practical significance to formulate optimization decisions to ensure economic benefits and navigation safety.

There are a total of eight turbines in XJB reservoir, four of which are located on the right-bank underground powerhouse, with the other four located on the left-bank power house at the dam toe, as shown in the Figure 3c. Besides, the left-bank turbines are about 900 m away from the entrance area of approach channel, with an interval of 36 m for each turbine. The right-bank turbines are about 1300 m away from the entrance area of approach channel, and each turbine is about 39 m apart.

In order to build the economic/navigation multi-objective model, the initial water level and the final water level of a day are both set as 379.5 m in this study. The average hourly inflow remains at 5000 m³/s. Moreover, the allowed navigation period of the XJB reservoir is from 08:00 to 18:00 h in a day, and a typical ship with a size of 85.0 m × 10.8 m × 2.0 m (length × width × draft depth) is selected for analysis. For the short-term optimization operation problem in this case study, the length of time in each period is set to one hour.

5.2. Parameter Settings and Simulation Working Conditions

5.2.1. Model Parameter Setting

According to the measured sections data of Jinsha River, two two-dimensional hydrodynamic models for the upstream and downstream reaches, respectively, are established by the software Mike21 FM to analyze the navigation capability. Moreover, the study regions of two models are selected, as shown in Figure 3d,e.

The width of the study region of the upstream reaches is 590 m and the length is 910 m. Besides, the total number of model grids is 1227, and the number of grid nodes is 690.

The size of the study area of the downstream model is $420\times 400$ m. Also, the total number of grid cells and the number of grid nodes is 1698 and 1008, respectively. The generated mesh of the upstream and downstream reservoir is shown in Figure 4.

Then, the numerical simulation results of the water level is compared with the measured data. In addition, the parameters of the model are modified by the result of the comparison until the simulated result matches the measured data within an acceptable range.

In order to explore the influence of the water level variation on navigation, the cross-section is set every 50 m in two models; finally, 18 cross-sections are set in the upstream and downstream models.

5.2.2. Simulation Working Conditions

According to the layout of turbines in Section 5.1, during the actual operation when all the turbines work normally, the symmetric flow will be generated in the downstream river reaches, because the discharge is uniformly distributed in the downstream reaches. In this condition, the upstream boundary line in the downstream model is AC, as shown in Figure 4, which means that the discharge is evenly distributed on the AC boundary.

However, if only some of the turbines work, then there will be asymmetric flow in the downstream river reaches because of the uneven discharge distribution. In order to explore the influence of asymmetric flow on the downstream navigation capacity, two typical working conditions are selected for research, as shown in

Table 2. Three simulation working conditions.

Condition	Discharge	Working Turbines	UBL 1	Simulation
1	Symmetric flow	All turbines	AC	Preliminary simulation and elaborate simulation
2	Asymmetric flow	Four left-bank turbines	AB	Asymmetric flow simulation
3	Asymmetric flow	Four right-bank turbines	BC	Asymmetric flow simulation

¹ UBL is the upstream boundary line of the downstream model.

.

Whether the symmetric flow or the asymmetric flow is determined by the outflow condition of the two sections AB and BC, as shown in Figure 4, which are on the two sides of the midpoint B of the boundary line on the downstream model. If there is outflow on both sections (AB and BC), a symmetric flow will be generated. If there is outflow on the only section (AB or BC), an asymmetric flow will be generated.

Thus, three simulations are performed in the case study, which are the preliminary simulation of symmetric flow, and the elaborate simulation of symmetric flow and asymmetric flow, respectively. Among them, Condition 1 is performed in the preliminary simulation and elaborate simulation, and Conditions 2 and 3 are completed in the asymmetric flow simulation.

5.3. Results

5.3.1. Numerical Simulation Results

(1) Preliminary simulation of symmetric flow

To briefly study the navigation capacity and determine the weighting factors (${\omega }_v$ and ${\omega }_l$) for the above two hydrodynamic models, four working conditions (two upstream conditions and two downstream conditions) are selected to simulate the influences of flow velocity and water level on ship navigation, as shown in

Table 3. The working conditions and results of preliminary simulation of symmetric flow.

Model	Condition	UBC 1	DBC 2	u (m/s) 4	v (m/s) 4	UBL
		Q (m3/s) 3	Z (m) 3
Upstream	P1	9880	380	0.08–0.12	0.04–0.05	——
	P2	9970	370	0.16–0.28	0.02–0.10	——
Downstream	P3	1200	265.8	0.05–0.35	0.025–0.125	AC
	P4	12,000	277.25	0.50–2.00	$\ge$0.3 (locally)	AC

¹ UBC is the upstream boundary condition; ² DBC is the downstream boundary condition; ³ Q is the inflow (I) and Z is the upstream water level (Z^up) in the upstream model; Q is the discharge flow (Qx), and Z is the downstream water level (Z^down) in the downstream model; ⁴ u is the longitudinal velocity and v is the transverse velocity.

. The value of discharge volume and water level are used to the upstream and downstream boundary conditions of the two models in the Table 3, respectively.

The upstream water level of Condition P1 is 380 m, which is the normal water storage level (September–December of this year, January and June of the following year). The highest flow occurs in September, and the average annual flow is 9880 m³/s. Besides, in Condition P2 the upstream water level is 370 m in July–September, the highest flow is in August, and the average annual flow is 9970 m³/s. Finally, Condition P3 and Condition P4 show the minimum and maximum navigable flow in the downstream reaches of the reservoir, respectively.

According to the model parameter settings listed in Section 5.2, the numerical simulation of the two hydrodynamic models is implemented separately. The research results are as follows:

(1)

The numerical simulation results of Condition P1 and Condition P2 are shown in Table 3, Figure 5, and Figure 6. It can be seen that the longitudinal velocity and transverse velocity of the upstream entrance area of approach channel satisfy the navigation requirements. Therefore, it can be inferred that the upstream navigation capacity ($N{C}^{up}$) is less affected by the flow velocity. In this study case, $N{C}^{up}$ is mainly related to the water level variation, and the formula for $N{C}^{up}$ is expressed as follows:

$$N{C}^{up}=\frac{1}{C{S}^{up}}\frac{1}{\left(Tn-1\right)}{\displaystyle \sum _{j=1}^{C{S}^{up}}{\displaystyle \sum _{t=2}^{Tn}(1-\left|\frac{Z_{t,j}^{down}-Z_{t-1,j}^{down}}{\Delta Z_{h,\lim \mathrm{t}}}\right|)}}$$

(22)

where: $C{S}^{up}$ = 12 and $\Delta Z_{h,\lim \mathrm{t}}$ = 1.5 m/h.

(2)

The transverse velocity and longitudinal velocity in Condition P3 and the longitudinal velocity in the Condition P4 are all satisfied by the navigation requirement, as shown in Figure 7 and Table 3. However, the local transverse velocity exceeds the limit value of 0.3 m/s in the Condition P4, which seriously hinders the navigation in the entrance area of the approach channel. As a result, the transverse velocity is the most important factor affecting the downstream navigation, which is considered in this case.

The downstream navigation should first consider the influence of the flow velocity and then consider the influence of the water level variation. Therefore, in the objective of the downstream navigation capacity ($N{C}^{down}$), the weight of the flow velocity accounts for a larger proportion (${\omega }_v$ = 0.95), and the water level variation accounts for a relatively small proportion (${\omega }_l$ = 0.05). Finally, the formula for $N{C}^{down}$ is expressed as follows:

$$N{C}^{down}=0.95\frac{1}{Tn}{\displaystyle \sum _{t=1}^{Tn}\left[1-\frac{1}{Dc}Dd(Q{x}_t,d_{ship})\right]}+0.05\frac{1}{C{S}^{down}}\frac{1}{\left(Tn-1\right)}{\displaystyle \sum _{j=1}^{C{S}^{down}}{\displaystyle \sum _{t=2}^{Tn}(1-\left|\frac{Z_{t,j}^{down}-Z_{t-1,j}^{down}}{\Delta Z_{h,\lim \mathrm{t}}}\right|)}}$$

(23)

where $C{S}^{down}$ = 22 and $Dd$ = 57 m.

(2) Elaborate simulation of symmetric flow

To explore the different values of $Dd$ and the cross-sectional water level under different discharge volume in the downstream model, 23 sets of numerical simulations with different discharge volumes ($Qx$) are implemented. The boundary conditions of simulations and the simulated results are shown in

Table 4. The results of numerical simulations in the downstream model.

No.	UBC	DBC	Dd	Cross-Sectional Water Level
	Qx	${\mathit{Z}}^{\mathit{d}\mathit{o}\mathit{w}\mathit{n}}$		No. 10	No. 11	No. 12	No. 13	No. 14	No. 15	No. 16	No. 17	No. 18
1	1200	265.8	0	265.80	265.80	265.80	265.80	265.80	265.80	265.80	265.80	265.80
2	1700	266.7	0	266.72	266.72	266.72	266.72	266.72	266.72	266.72	266.72	266.72
3	2200	267.6	0	267.56	267.56	267.56	267.56	267.56	267.56	267.56	267.56	267.56
4	2700	268.3	0	268.29	268.29	268.29	268.29	268.29	268.29	268.29	268.29	268.29
5	3200	269.0	0	269.00	269.00	269.00	269.00	269.00	269.00	269.00	269.00	269.00
6	3700	269.6	0	269.63	269.63	269.63	269.63	269.63	269.63	269.63	269.63	269.63
7	4200	270.2	25	270.24	270.24	270.24	270.24	270.24	270.24	270.24	270.24	270.24
8	4700	270.8	28	270.82	270.82	270.82	270.82	270.82	270.82	270.82	270.82	270.82
9	5200	271.3	31	271.35	271.35	271.35	271.35	271.35	271.35	271.35	271.35	271.35
10	5700	271.8	34	271.84	271.84	271.84	271.84	271.84	271.84	271.84	271.84	271.84
11	6200	272.3	36	272.32	272.32	272.32	272.32	272.32	272.31	272.31	272.31	272.31
12	6700	272.8	37	272.77	272.77	272.77	272.77	272.77	272.77	272.76	272.76	272.75
13	7200	273.2	38	273.21	273.21	273.21	273.21	273.21	273.21	273.20	273.20	273.19
14	7700	273.6	40	273.64	273.64	273.64	273.64	273.64	273.64	273.63	273.63	273.62
15	8200	274.1	41	274.07	274.07	274.07	274.07	274.07	274.07	274.06	274.06	274.05
16	8700	274.5	42	274.49	274.49	274.49	274.49	274.49	274.49	274.48	274.48	274.47
17	9200	274.9	43	274.92	274.92	274.92	274.92	274.92	274.92	274.91	274.91	274.90
18	9700	275.3	43.5	275.34	275.34	275.34	275.34	275.34	275.34	275.33	275.32	275.31
19	10,200	275.8	43.6	275.77	275.77	275.77	275.77	275.77	275.77	275.76	275.75	275.74
20	10,700	276.2	43.7	276.19	276.19	276.19	276.19	276.19	276.19	276.18	276.17	276.16
21	11,200	276.6	43.8	276.60	276.60	276.60	276.60	276.60	276.59	276.59	276.58	276.57
22	11,700	277.0	43.9	277.02	277.02	277.02	277.02	277.02	277.01	277.01	277.00	276.98
23	12,000	277.3	44.5	277.27	277.27	277.27	277.27	277.27	277.26	277.26	277.25	277.23

. According to the water level–flow relationship of XJB reservoir, the discharge volume ($Qx$) profile is applied as the upstream boundary condition, and the downstream water level ($Z^{down}$) profile is used as the downstream boundary condition in the downstream model.

Taking the discharge volume of 5700 m³/s as an example, $NC{v}^{down}$ = 0.403 is obtained through the calculation, as shown in Figure 8.

In order to calculate the value of cross-sectional water level under different conditions in the upstream model, 21 sets of numerical simulations are completed, and the results are shown in

Table 5. The results of numerical simulations in the upstream model.

No.	UBC	DBC	Cross-Sectional Water Level
	I	${\mathit{Z}}^{\mathit{u}\mathit{p}}$	No. 1	No. 2	No. 3	No. 4	No. 5	No. 6	No. 7	No. 8	No. 9
1	5000	380.0	380.03	380.03	380.02	380.02	380.02	380.01	380.01	380.01	380.00
2	5000	379.5	379.50	379.50	379.50	379.50	379.50	379.50	379.50	379.50	379.50
3	5000	379.0	379.00	379.00	379.00	379.00	379.00	379.00	379.00	379.00	379.00
4	5000	378.5	378.50	378.50	378.50	378.50	378.50	378.50	378.50	378.50	378.50
5	5000	378.0	378.00	378.00	378.00	378.00	378.00	378.00	378.00	378.00	378.00
6	5000	377.5	377.50	377.50	377.50	377.50	377.50	377.50	377.50	377.50	377.50
7	5000	377.0	377.00	377.00	377.00	377.00	377.00	377.00	377.00	377.00	377.00
8	5000	376.5	376.50	376.50	376.50	376.50	376.50	376.50	376.50	376.50	376.50
9	5000	376.0	376.00	376.00	376.00	376.00	376.00	376.00	376.00	376.00	376.00
10	5000	375.5	375.50	375.50	375.50	375.50	375.50	375.50	375.50	375.50	375.50
11	5000	375.0	375.00	375.00	375.00	375.00	375.00	375.00	375.00	375.00	375.00
12	5000	374.5	374.50	374.50	374.50	374.50	374.50	374.50	374.50	374.50	374.50
13	5000	374.0	374.00	374.00	374.00	374.00	374.00	374.00	374.00	374.00	374.00
14	5000	373.5	373.50	373.50	373.50	373.50	373.50	373.50	373.50	373.50	373.50
15	5000	373.0	373.00	373.00	373.00	373.00	373.00	373.00	373.00	373.00	373.00
16	5000	372.5	372.50	372.50	372.50	372.50	372.50	372.50	372.50	372.50	372.50
17	5000	372.0	372.00	372.00	372.00	372.00	372.00	372.00	372.00	372.00	372.00
18	5000	371.5	371.50	371.50	371.50	371.50	371.50	371.50	371.50	371.50	371.50
19	5000	371.0	371.00	371.00	371.00	371.00	371.00	371.00	371.00	371.00	371.00
20	5000	370.5	370.50	370.50	370.50	370.50	370.50	370.50	370.50	370.50	370.50
21	5000	370.0	370.00	370.00	370.00	370.00	370.00	370.00	370.00	370.00	370.00

. According to the average hourly inflow (I) set in Section 5.1, the inflow profile is used as the upstream boundary condition, and the upstream water level ($Z^{up}$) profile is applied as the downstream boundary condition in the upstream model. Therefore, the water levels of the downstream boundary condition are selected between the dead water level (370 m) and the normal height (380 m), with an interval of 0.5 m.

(3) Asymmetric flow simulation

To study the influence of asymmetric flow on the $NC{v}^{down}$ of XJB reservoir, Conditions 2 and 3 in Table 2 are simulated by the same downstream model with different boundary lines.

Twenty-three sets of numerical simulations are carried out, and their boundary conditions are the same as those in the elaborate simulation of symmetric flow, as shown in Table 4. The calculation results are shown in Figure 9.

Through the simulation results, the following conclusions can be drawn:

(1)

For downstream navigation of XJB reservoir, when the discharge volume (less than 1700 m³/s) is relatively small, the downstream navigation capacity considering the flow velocity ($NC{v}^{down}$) of the three working conditions is the same ($NC{v}^{down}$ = 1);

(2)

When the discharge volume is large (greater than 1700 m³/s), the $NC{v}^{down}$ is the poorest in Condition 1 (only left-bank turbines work), better in Condition 2 (all turbines work), and the best in Condition 3 (only right-bank turbines work).

Due to the fact that the left-bank turbines are close to the entrance area of approach channel, if only the left-bank turbines are working, then a larger area where the transverse velocity exceeds the limit value will be generated, and there will be a greater dangerous navigation width ($Dd$). In contrast, the right-bank turbines are far away from the entrance area and have less impact on the navigation of the entrance area. Therefore, the study and results of asymmetric flow can provide useful suggestions for the unit commitment problem [45]. For example, when less power generation is required by XJB reservoir, it is preferable to only operate the right-bank turbines in order to satisfy the navigation demands.

The length of time in each period is set to 1 h in this paper, so the symmetric flow is selected as the study object in the short-term optimal operation problem. Besides, the study results of asymmetric flow can provide guidance for future ultra-short term optimal operation problems or inner-plant economical operations [46].

5.3.2. Multi-Objective Model Results

Based on the simulation results of symmetric flow in the two models (Table 4 and Table 5) and the determination of weighting factors (${\omega }_v$ and ${\omega }_l$), the objective functions of $N{C}^{up}$ and $N{C}^{down}$ are expressed by Formulas (22) and (23), respectively. Then, the multi-objective optimal operation model is calculated according to the process of Figure 2. The algorithm parameters are set as follows: the population size is set to 200, the size of the archive is set to 100, the number of iterations is set to 20,000 generations, the mutation rate is set to 0.04, and the crossover rate is set to 0.85.

The Pareto optimal front of three objectives is acquired as given in Figure 10a; the optimal results of any two of the objectives are also expressed in Figure 10b–d. All the operation schemes results of the reservoir are shown in

Table 6. Operation schemes of reservoir (arranged in ascending order of $N{C}^{down}$).

Scheme	$\mathit{E}$ (104 kWh)	$\mathit{N}{\mathit{C}}^{\mathit{d}\mathit{o}\mathit{w}\mathit{n}}$	$\mathit{N}{\mathit{C}}^{\mathit{u}\mathit{p}}$	Scheme	$\mathit{E}$ (104 kWh)	$\mathit{N}{\mathit{C}}^{\mathit{d}\mathit{o}\mathit{w}\mathit{n}}$	$\mathit{N}{\mathit{C}}^{\mathit{u}\mathit{p}}$
1	10,942.07743	0.516103	0.990381	40	10,934.87	0.775069	0.982459
2	10,942.65972	0.522835	0.988732	41	10,938.68	0.776607	0.981775
3	10,942.33743	0.538384	0.988012	42	10,933.95	0.778997	0.982459
4	10,942.11842	0.550021	0.989657	43	10,937.91	0.783004	0.978221
5	10,941.27857	0.553263	0.992386	44	10,934.97	0.790506	0.982092
…	…	…	…	…	…	…	…
35	10,939.44119	0.742237	0.981028	96	10,922.86	0.979242	0.968444
36	10,939.48026	0.745708	0.980875	97	10,925.08	0.987635	0.96826
37	10,938.9987	0.750805	0.979485	98	10,924.49	0.989242	0.968253
38	10,938.77714	0.758043	0.979526	99	10,923.75	0.989259	0.968254
39	10,938.71673	0.763488	0.978766	100	10,925.44	0.99352	0.968062

. The three conclusions drawn from the Figure 10 are as follows:

(1)

There is an obvious inverse relationship between the total power generation ($E$) and the downstream navigation capacity ($N{C}^{down}$) from the results of Figure 10b. The larger the total power generation ($E$), the smaller the downstream navigation capacity ($N{C}^{down}$) in a day. The minimum value of $N{C}^{down}$ is 0.52, and the maximum value is 0.99 with a growth of 90.38%, which varies greatly. At the same time, there is a drop of $E$ from the maximum value at 10,942.07 × 10⁴ kWh to the minimum value at 10,925.44 × 10⁴ kWh.

(2)

As shown in the Figure 10c, there is a certain inverse trend between the upstream navigation capacity ($N{C}^{up}$) and the downstream navigation capacity ($N{C}^{down}$). When the upstream navigation capacity increases, the downstream navigation capacity declines. The minimum value of $N{C}^{up}$ in all the schemes is 0.968, and the maximum value is 0.992 with a smaller growth compared to the change in $N{C}^{down}$.

(3)

Finally, it can be seen that the relationship between the total power generation ($E$) and the upstream navigation capacity ($N{C}^{up}$) is not obvious shown in the Figure 10d. There is little interaction between these two elements.

5.4. Discussion

To further understand the relationship between three objectives and explain the cause of this phenomenon, typical schemes are selected to be discussed and analyzed in depth. Because $N{C}^{up}$ and $E$ have no obvious relationship and the change of $N{C}^{up}$ in all schemes is relatively small, the three schemes which mainly consider the relationship between $N{C}^{down}$ and $E$ are selected. The schemes consist of the maximum of the total power generation (Scheme 1), the extreme point of the maximum of the downstream navigation capacity (Scheme 100), and the relative intermediate point (Scheme 37) compared to other schemes. The comparison of the objective values of the three schemes is expressed in Figure 11a, and the specific power generation and navigation processes of the three schemes are shown in Figure 11b–d, respectively.

The common feature of the three schemes is that the water level drops from 00:00 h to 0:00 h (the start of navigation), and is in a state of gradual drainage before navigation. Then the reservoir begins to store the water and the water level grows between 08:00 and 18:00 h (the end of navigation), thus decreasing the amount of the discharge volume in the navigable period in order to establish a favorable navigation environment. After the navigation ends, the water level drops to the end water level in the following hours, which will maximize the amount of power generation after the navigation process.

In Scheme 1, the water level falls from 00:00 to 08:00 h to 379.37 m slowly, with a moderate drop, followed with a slight rise during the navigation time. However, during the navigation period, the average discharge volume of the reservoir is maintained at 4602 m³/s, and the average output is maintained at 4213 MW, resulting in a large amount of electricity during this period. In the meantime, the value of $NC{v}^{down}$ per hour is maintained between 0.48 and 0.52, resulting in a large transverse velocity in the entrance area of approach channel, which is not conducive to navigation according to the numerical simulation results in Section 5.3.1. This scheme significantly rises the total power generation, but reducing the downstream navigation capacity in contrast.

In Scheme 100, after the water level sharply declines from 00:00 to 08:00 h, the reservoir strengthens the water storage, and the water level increases rapidly. During the navigation time, the average output and discharge volume are maintained at 3430 MW and 3708 m³/s, respectively. Consequently, the value of $NC{v}^{down}$ per hour is also the highest with an average value of 1, and it builds beneficial condition for navigations while decreasing the economic benefits of the reservoir.

Scheme 37 represents a compromise that considers both economics and navigation. It can provide short-term power generation operation advice for reservoir managers, and it can also guide ship operators to safely pass through dangerous areas.

6. Conclusions

In the past, most navigation capacity evaluation methods neglected the influence of flow velocity on navigation of ships. Based on the two-dimensional hydrodynamic numerical simulation, a navigation capability evaluation method that considers the flow velocity and water level comprehensively is proposed in this paper. In addition, a short-term multi-objective optimal operation model considering the upstream and downstream navigation and power generation of the reservoir is established and solved by SPEA2. To verify the reasonableness and practicability of the method and model, an example of short-term optimal operation of the Xiangjiaba reservoir in China is studied. Based on the case study, this paper analyzes in depth the impact of various operation processes of reservoir on economic and navigable benefits. The following conclusions can be drawn:

(1)

The proposed NCEM to evaluate the navigation capability is reasonable and effective, and can comprehensively analyze the influence of flow velocity and water level variation on navigation accurately.

(2)

The proposed multi-objective model can obtain a favorable Pareto frontier and explore the relationship between objectives. In the case study of the XJB reservoir, there is an obvious inverse relationship between power generation and the downstream navigation capacity. Also, the relationship between downstream navigation capacity and upstream navigation capability is inverse. However, there is little interaction between power generation and upstream navigation capability.

(3)

The results illustrate that the method and model are reasonable and effective, and also indicate that they can provide a series of favorable optimal operation schemes for the reservoir to obtain economic and navigational benefits.