1 Introduction

Hydropower-generated energy is an indispensable and sustainable energy source. By harnessing the power of flowing water, this renewable energy option offers numerous advantages, especially regarding the current environmental changes happening on the planet. Hydroelectric power plants (HPPs) provide a responsible and eco-friendly solution for generating electricity while minimizing environmental impact. However, one of the greatest challenges lies in effectively modeling the operations of these power plants, considering both power generation and their associated impacts. Extensive studies have been conducted on the environmental effects of HPPs [1], and methodologies have been developed to mitigate these impacts [2, 3].

Many approaches consist of optimizing an energy dispatch model, Operation Multi-Reservoir System (OMRS) [4, 5]. In the context of multi-reservoir systems, the optimization problem presents inherent challenges, such as a non-linear objective function and a vast search space. Not rarely, multi-objective models are also considered. Methods for solving OMRS optimization problems vary from classical deterministic methods to stochastic evolutionary meta-heuristics. Different traditional methods have been employed in the optimization of HPPs’ operations, such as Linear Programming (LP) [6], Mixed Integer Linear Programming (MILP) [7, 8] and Non-Linear Programming (NLP) [9, 10].

Bio-inspired meta-heuristics have also been employed as a method for solving OMRS problems. In [11, 12], the authors make use of evolutionary algorithms to optimize energy generation in hydropower plants. The work by [13] uses Particle Swarm Optimization (PSO) to find the optimal benefit-cost ratio in Hydro Power Plants design, PSO was also applied to the optimization of daily scheduling in HPPs reservoirs [14]. Multi-objective approaches have gained prominence in recent years. As examples we can highlight, a combination approach for downstream plants to solve scheduling information asymmetry problem in electricity markets was realized by [15], an optimization of long-term operation of the large Karun in Iran was carried out by [16], and a study on the regulation potential of the upper Yellow River reservoirs in China discussed by [17].

Even though evolutionary algorithms have presented a certain degree of success in the electric dispatch problem, there are specific weaknesses associated with these approaches, especially related to their non-deterministic nature. Such methods are not guaranteed to find the global optima and at times, they get trapped in local optima due to premature convergence [18, 19]. Besides the issue of local optima, it is common that problems present a multimodal structure – the case where different parts of the Pareto Frontier are scattered across the decision space. Recently, many studies have been conducted on techniques to address multimodal multi-objective problems [20, 21], as well as theoretical studies which aim to provide a better understanding of these search spaces [22].

In this work, we investigated the impact of including the Special-Crowd-Distance operation in the MESH algorithm and verified the ability of this algorithm to solve a problem with three conflicting objectives, which had not yet been studied. The paper has been organized as follows: Sect. 2 describes the model of the operation multi-reservoir system with three conflicting objectives. Section 3 details the MESH algorithm proposed using the Special-Crowd-Distance operator by [21]. The simulation, results and comparative analysis of algorithms are shown in Sect. 4. Conclusions regarding the achieved performance of Special-Crowd-Distance boosted MESH algorithm are presented in Sect. 5.

2 Mathematical Model for the Operation of Cascade Hydro-Power Plants

Hydropower plants typically operate using one of two approaches: the Joint Control System and the Optimized Control System. The Joint Control System, where the amount of water flow is divided equally among every available turbine generator, and the Optimized Control System, where the amount of water flow in each turbine generator is decided to maximize the power generation. The model described in the current section assumes an Optimized Control System, where the decision space of the optimization model is the water flow in each of the generating units of a hydro-power plant.

Traditionally, power generation in hydropower plants has been optimized based solely on maximizing energy output. However, as sustainability and environmental concerns have gained prominence, the need for a more comprehensive approach has become evident. A multi-objective model takes into account not only the energy generation potential but also other factors such as ecological flow requirements, downstream water quality, and reservoir management. The 2-objective model proposed in [23] and [12] aims to maximize the power generated by the cascade system and the volume of water in each of the HPPs reservoirs. In this work such a model is improved by an objective function, proposed in [24], which accounts for a better human control of the optimized flow in each turbine generator. As a result, the model is classified as a 3-objective optimization model.

Table 1. Parameter table of the cascade HPP model
Full size table

In such electric dispatch model, the power generated at each timestamp for generating unit (j) at HPP (u), can be calculated with Eq. 1 (See Table 1 for a description of the variables).

$$\begin{aligned} \begin{aligned} ph_{u,j,t} = g \cdot [ \rho _{0} + \rho _1 hl_{u,j,t} + \rho _2 q_{u,j,t} + \rho _3 hl_{u,j,t} q_{u,j,t} \\ + \rho _4 hl_{u,j,t}^2 + \rho _5 q_{u,j,t}^2] \cdot [ hb_{u,j,t} - \varDelta _{Hu,j,t}] \cdot q_{u,j,t}. \end{aligned} \end{aligned}$$
(1)

In which, \(\rho _1, \rho _2, \rho _3, \rho _4,\rho _5\) are coefficients of the power generation equation estimated through ordinary least squares method from historical data, and \(\mathbf {q_t} = [q_{1,1,t}, q_{1,2,t}, ..., q_{U,J_u, t}]\) is the water flow in each turbine generator, i.e. the decision space.

Moreover, the values of essential variables in the calculation of Eq. 1 can be obtained with the following equations

$$\begin{aligned} hl_{u,j,t} = hb_{u,j,t} - \varDelta _{Hu,j,t}, \end{aligned}$$
$$\begin{aligned} hb_{u,j,t} = fcm_{u,t} - fcj_{u,t}, \end{aligned}$$
$$\begin{aligned} fcm_{u,t} = a_0 + a_1 \cdot \varPsi _{u, t} + a_2 \cdot \varPsi _{u, t}^2 + a_3 \cdot \varPsi _{u, t}^3 + a_4 \cdot \varPsi _{u, t}^4, \end{aligned}$$
(2)
$$\begin{aligned} \begin{aligned} fcj_{u,t} = b_0 + b_1 \cdot (Qt_{u,t} + Qv_{u,t}) + b_2 \cdot (Qt_{u,t} + \\ Qv_{u,t})^2 + b_3 \cdot (Qt_{u,t} + Qv_{u,t})^3 + b_4 \cdot (Qt_{u,t} + Qv_{u,t})^4, \end{aligned} \end{aligned}$$
(3)

as shown in Table 1, \(fcm_{u,t}\) (Eq. 2) and \(fcj_{u,t}\) (Eq. 3) are, respectively, the height of the HPP (u) downstream and upstream at time (t), and coefficients \(a_0, a_1, ..., a_4\) and \(b_0, b_1, ..., b_4\) are estimated with the ordinary least squares method from historical data.

In a cascade HPPs context, the water balance in the reservoir of each power plant plays a fundamental role, since in such a scenario, the reservoir of a single HPP is influenced by the flow rates of all the previous HPPs in the system. The water balance of each HPP in the system is, also, which accounts for the dynamical nature of the optimization problem, since the volume of water in the reservoir is time dependent. In the model proposed in [12] the volume of water in HPP (u) at time t is defined by Eq. 4

$$\begin{aligned} \begin{aligned} \varPsi _{u,t} = \varPsi _{u,t - 1} + Qa_{u, t-1} + Qt_{w,td} + Qv_{w,td} - \\ Qt_{u,t-1} - Qv_{u,t-1} - E_{u,t-1} \cdot A_{u,t-1}. \end{aligned} \end{aligned}$$
(4)

In which, \(\varPsi _{u,t}\) is the reservoir volume of HPP u at time t, (u) and (w) are distinct HPPs in the system and (td) is the time necessary for the flow to go from (w) into (u). Additionally, Qa denotes the affluent flow, Qt represents the turbine flow, Qv represents the total flow rate, E corresponds to liquid evaporation, and A represents the area occupied by water in the reservoir (See Table 1).

The first objective function, denoted as \(F_1\) and proposed in [23], refers to the energy efficiency of the model and aims to maximize the power generated \(ph_{jt}\) per unit water flow. It is formulated as follows:

$$\begin{aligned} \text {Maximize} \; F_1(\mathbf {q_t}) = \frac{1}{U} \sum _{u=1}^U\left( \frac{\sum _{j=1}^{J_u} ph_{uj,t}}{\sum _{j=1}^{J_u} q_{uj,t}}\right) . \end{aligned}$$

Here, U represents the number of cascaded power plants, and \(J_u\) denotes the number of generating units in power plant u.

The second objective function, \(F_2\), proposed in [12], aims to maximize the water levels in the reservoir. Maintaining high water volume levels enhances the robustness of the hydropower plant during dry periods. Moreover, by increasing the water volume, the value of the upstream level (cota montante) is also elevated, thus increasing the power output. The function is given by:

$$\begin{aligned} \text {Maximize} \; F_2(\mathbf {q_t}) = \frac{1}{U} \sum _{u=1}^U\left( \frac{V_{u,t}}{V_u^{max}}\right) . \end{aligned}$$

where \(V_u^{max}\) represents the maximum reservoir volume.

The third objective function, denoted as \(F_3\), is adapted for cascade hydro-power plants from [24]. It minimizes the Euclidean distance between the optimized flow rates of the generating units and the traditional flow rates in the joint control system, where all generating units have the same flow rate (as is typically operated). This objective aims to increase the confidence of the hydro-power plant operators in the optimized flow rates by seeking those that are closest to the coordinated control flow rates while maximizing energy efficiency and reservoir water volume. The function is formulated as:

$$\begin{aligned} \text {Minimize}\; F_3(\mathbf {q_t}) = \frac{1}{U} \sum _{u=1}^U \left( \sum _{j=1}^{J_u} (q_{uj,t} - q_{u,cc})^2\right) . \end{aligned}$$

Here, \(q_{u, cc}\) represents the flow rate of the generating units in power plant u in the coordinated control model.

In addition to the objective functions, the decision space is subject to certain constraints [12]. The first constraint, described in Eq. 5, ensures that the energy generated by the hydro-power plant matches the demand \(Dm_{u,t}\) with a specified error tolerance \(\epsilon \) (approximately 0.5%).

$$\begin{aligned} Dm_{u,t} \cdot (1 - \epsilon ) \le \sum _{j=1}^{J_u} ph_{uj,t} \le Dm_{u,t} \cdot (1 + \epsilon ). \end{aligned}$$
(5)

Furthermore, the water flow in each generating unit must comply with the minimum flow rate \(q_{uj}^{min}\) and maximum flow rate \(q_{uj}^{max}\) limits of power plant u. The same applies to the water volume in the reservoir, which should fall within the minimum volume \(V_u^{min}\) and maximum volume \(V_u^{max}\) limits, according to:

$$\begin{aligned} q_{uj}^{min} \le q_{uj,t} \le q_{uj}^{max}, \end{aligned}$$
$$\begin{aligned} V_u^{min} \le V_{u,t} \le V_u^{max}. \end{aligned}$$

The defluent flows of the hydro-power plants also need to respect the interval \([Qd_{u}^{min}, Qd_{u}^{max}]\) to prevent flooding in the regions surrounding the plants, as:

$$\begin{aligned} Qd_{u}^{min} \le Qv_{u,t} + \sum _{j=1}^{J_u} q_{uj,t} \le Qd_{u}^{max}. \end{aligned}$$

A maximum limit \(Qv_u^{max}\) is imposed on the spill flow rate as well:

$$\begin{aligned} Qv_{u, t} \le Qv_u^{max}. \end{aligned}$$

Finally, the power generated by generating unit j must fall within the minimum \(ph_{uj}^{min}\) and maximum \(ph_{uj}^{max}\) limits for that unit in power plant u,

$$\begin{aligned} ph_{uj}^{min} \cdot Z_{uj,t} \le ph_{uj,t} \le ph_{uj}^{max} \cdot Z_{ujt} . \end{aligned}$$

In which, \(Z_{uj,t} \in {0, 1}\) is a binary variable indicating whether generating unit j of power plant u is active (\(Z_{uj,t} = 1\)) or inactive (\(Z_{uj,t} = 0\)) at time t.

3 Multi-objective Evolutionary Swarm Hybridization Algorithm

Evolutionary algorithms are a class of computational techniques inspired by the process of biological evolution. They are widely used in optimization problems, particularly in cases where traditional mathematical approaches are difficult to apply. Evolutionary algorithms mimic the principles of natural selection, reproduction, and mutation to iteratively search for optimal or near-optimal solutions. When it comes to multi-objective problems, where multiple conflicting objectives need to be simultaneously optimized, evolutionary algorithms offer a powerful approach. One of such is the Multi-objective Evolutionary Swarm Hybridization algorithm (MESH) proposed in [12].

The algorithm is a multi-objective version of the previously proposed C-DEEPSO algorithm [25], which implements hybridization techniques by combining operators from Differential Evolution [26] and Particle Swarm Optimization [27]. As in many optimization approaches through evolutionary algorithms, MESH aims to find the optimal solution by keeping a population of solutions each iteration, which is updated according to the evolutionary operators. The swarm optimization methodology, followed by MESH, finds the optimal solution by updating every particle with a velocity vector that captures information from the entire population of where the optimal solution is. This procedure, for a single particle \(X_i^t\) is illustrated in Eqs. 6 and 7.

$$\begin{aligned} X^{t + 1}_{i} = X^{t}_i + V^{t}_i, \end{aligned}$$
(6)
$$\begin{aligned} V^{t}_i = w_I^* \cdot V^{t - 1}_i + w_A^* \cdot (X^{t}_{sn} - X^{t}_i) + w_C^* \cdot C \cdot (X_{gb}^* - X^{t}_i). \end{aligned}$$
(7)

Similar to the procedure carried out by Particle Swarm Optimization (PSO), the velocity vector \(V^{t}_i\) is constructed by the sum of three components, named as inertia, cognitive and social components. The main difference between the procedure described in Eq. 7 and PSO’s resides in the cognitive component, where the term \(X^{t}_{sn}\) is obtained by applying the recombination mechanism from Differential Evolution (DE).

The superscript * on several parameters in Eq. 7, indicates that such parameter suffers a stochastic mutation process, given by Eq. 8.

$$\begin{aligned} w^* = w + \tau \cdot N(0,1). \end{aligned}$$
(8)

In which \(\tau \) is one of the algorithm’s hyperparameters and N(0, 1) is a sample from a Standard Gaussian Distribution. To account for the multi-objective nature of the problems, the algorithm stores a subset of the best particles found so far. Each iteration such subset of size N, located in a separate memory file, is updated through the crowd-distance operator from NSGA-II [28].

The crowd-distance operator from NSGA-II [28] begins by ordering each of the particles \(X_i\) in the population into dominance ranks as follows: in the first rank \(F_1\), there are all the particles not dominated by any other particle in the population P; in the second rank \(F_2\), there are those not dominated by any other particle in the population \(P - F_1\); in \(F_3\), there are those not dominated by any other particle in the population \(P - F_1 - F_2\), and so on.

After that, for each particle \(X_i\), its crowd-distance \(I_{i}\) is calculated based on the average distances between the two adjacent particles in all dimensions of the objective space [28, 29].

$$\begin{aligned} \begin{aligned} I_{i}^j = \frac{X^{i + 1}_j - X^{i - 1}_j}{f_j^{\text {max}} - f_j^{\text {min}}}, I_{i} = \sum _{j = 1}^n I_i^j. \end{aligned} \end{aligned}$$

In which, \(X^{i + 1}_j\) and \(X^{i - 1}_j\) are the values of the j-th objective function of the two adjacent particles to \(X^i\) when the particles are sorted based on their values for the j-th objective function. \(f_j^{\text {max}}\) and \(f_j^{\text {min}}\) are the maximum and minimum values for such function, respectively.

Then, the highest ranks (\(F_1, F_2, ...\)) are directly added to the memory file until the entire space of N particles is occupied. If only part of a rank can be inserted into memory, we sort the particles in descending order according to their crowd-distance, and these particles are added in order until the memory is completely full.

MESH also admits different configurations for how it operates, those are:

  • E, the choice of operator for the PSO swarm guide.

    • (E1) Select the nearest particle in memory.

    • (E2) Select the nearest particle in the next upper boundary to the current particle’s boundary.

  • V, from which group the particles that compose the differential evolution operator come from.

    • (V1) Population.

    • (V2) Memory.

    • (V3) Population + Memory.

  • D, which of the strategies for the differential evolution operator, will be used.

    • (D1) DE/Rand/1/Bin.

    • (D2) DE/Rand/2/Bin.

    • (D3) DE/Best/1/Bin.

    • (D4) DE/Current-to-best/1/Bin.

    • (D5) DE/Current-to-rand/1/Bin.

Therefore the total amount of possible different ways MESH can operate is \(2 \times 3 \times 5 = 60\). In this work, we focus on the E2V2D1 configuration, which, from empirical tests, appears to be the most stable one for general problem-solving.

3.1 Special Crowd Distance

In addition to the original configurations of the Multi-objective Evolutionary Swarm Hybridization (MESH) algorithm, this work implemented an additional parameter related to the crowd distance operator.

Unlike the original crowd distance operator defined by the NSGA-II algorithm, the additional configuration allows for the use of the special crowd distance operator defined in [21].

Many optimization problems exhibit the characteristic of being multimodal, meaning they have distinct local maxima in different regions of the decision space. In this context, it is important to consider diversity not only in the objective space (as in the original NSGA-II operator) but also in the decision space. To calculate the crowd distance in the decision space, we use the following expression proposed by [21]:

$$\begin{aligned} CD_{i,x} = \sum _{j = 1}^{k} (k - j + 1) d_{i, j}. \end{aligned}$$

In which, \(CD_{i,x}\) is the crowd distance in the decision space, i is the index of a particle in the population, j is the index of a particle in the neighborhood of the particle with index i, \(d_{i,j}\) is the Euclidean distance from the i-th particle in the entire population to the j-th particle in the neighborhood, and k is the neighborhood size.

The crowd distance in the objective space of the i-th particle, denoted as \(CD_{i, f}\), is calculated as in the original NSGA-II operator with a slight modification: the crowd distance is not calculated considering only one Pareto rank but the entire population [21].

Thus, the special crowd distance of the i-th particle, denoted as \(SCD_i\), is obtained, according [21] :

$$\begin{aligned} \begin{aligned} SCD_{i} = {\left\{ \begin{array}{ll} max(CD_{i,x}, CD_{i,f}), & \text {if } CD_{i,x} > CD_{avg, x} \\ \text { or } CD_{i,f} > CD_{avg, f}\\ min(CD_{i,x}, CD_{i,f}), & \text {otherwise} \end{array}\right. } \end{aligned} \end{aligned}$$

where \(CD_{avg, x}\) and \(CD_{avg, f}\) are the averages of the crowd distances in the decision space and the objective space, respectively. In the special crowd distance operator, the particles in the last rank to enter the memory are chosen according to the ordering of \(SCD_i\).

To incorporate the special crowd distance [21] into the MESH algorithm, a new configuration for MESH is added.

  • (C1) Traditional crowd distance operator from NSGA-II algorithm.

  • (C2) Special Crowd Distance from [21].

And the total number of possible configurations for MESH, with the crowd distance new parameter, is \(2 \times 3 \times 5 \times 2 = 120\). A previous study was carried out in [12]. Then, in this work the E2V2D1–MESH configuration was used for tested Special Crowd Distance operation.

4 Experiments and Results

The experiments conducted in this study compare the performance of the Special-Crowd-Distance boosted MESH algorithm (E2V2D1C2 version) against three other well-established metaheuristics: standard MESH (E2V2D1C1 version) [12], NSGA-II [28], NSGA-III [30, 31], and MOEA/D [32]. All three algorithms were implemented in Python version 3.7.12 using the pymoo library [33]. The MESH algorithm was evaluated in two distinct versions: E2V2D1C1, which uses the NSGA-II crowd distance operator, and E2V2D1C2, which incorporates a special crowd distance operator based on [21].

For all four algorithms, the termination criterion was set to 15.000 objective function evaluations and a population size of 100 particles. The algorithms are tested on their performance in solving a real-world problem concerning electric dispatch and were executed 30 times, as described in Sect. 2.

The main criterion employed for comparing the performance of the algorithms in the experiments is hypervolume, a widely used metric in multi-objective optimization. To assess the statistical significance of the observed results in hypervolume values among different algorithms, the Mann-Whitney U test [34] is employed. The Mann-Whitney U test is used to determine whether two independent samples originate from the same distribution or if there is a significant difference between them. Such a test was chosen due to its non-parametric nature, as the hypervolume samples do not appear to be drawn from a Gaussian distribution. By utilizing this test, we can confidently evaluate the performance disparities between the MESH algorithm versions and the other meta-heuristics under consideration.

Moreover, the statistical test was applied between the algorithm with the highest hypervolume results and each of the other algorithms. By conducting these pairwise comparisons, we can determine which algorithm exhibits the best performance in terms of hypervolume. The intention behind this methodology is to appropriately declare the algorithm that outperforms the others based on the statistical significance of the results.

The proposed model for the cascade hydro-power plant is dynamic, and a strategy is employed to reevaluate the Pareto Front as it changes over time. For each hour, 30 algorithm runs are performed, combining the solutions found into a combined Pareto Front. The Fig. 1 shows the combined Pareto Front obtained by each algorithm in a sample of hours of a day. In the context of Combined Pareto Frontiers, MOEA/D did not perform satisfactorily, its results for the hypervolume metric are consistently lower than those of the other algorithms in all hours of operation of the Hydroelectric Plants. The best results for all hours were between the MESH versions (E2V2D1C1 and E2V2D1C2), with a certain tendency towards the E2V2D1C2 configuration that we are proposing in this work, which obtained the highest results for hypervolume in most HPPs operating hours.

In general, 4 of the 5 algorithms were able to find consistent Pareto Frontiers and proved to be efficient when applied to the optimization of the electrical dispatch model in OMRS scheme. Except for MOEA/D, all algorithms appear to have good performance in finding the Pareto frontier for all hours of operation of the hydroelectric plant. Finally, it is notable that the MESH versions better cover the range of the solution set, at all hours of operation. From this set of solutions, inside the combined Pareto Front, a single solution must be selected. The values of the decision space variables in this selected solution will determine the future values of the dynamic variables. In this work, the selected solution is the one closest to the point in the middle of Pareto Front.

Table 2 contains the results for the average and standard deviation of hypervolume in each of the tested algorithms. E2V2D1C2 configuration of MESH presented the highest average values for most operating hours. However, its results were very close to the results of the E2V2D1C1 configuration, in a way that no statistically significant difference was detected in 17 out of 24 h. Moreover, both MESH configurations achieved superior results to the well-established algorithms in the literature. In addition to the good hypervolume results for both MESH versions, in Table 2, it is noticeable that the standard deviation results for both configurations are consistently more controlled for all operating hours compared to the other three algorithms evaluated in this experiment.

Fig. 1.
figure 1

Combined Pareto Fronts

Full size image
Table 2. Average and standard deviation (in bracket) results for the hypervolumes of each algorithm in daily electric dispatch per hour (H). The symbol ‘=’ means there are no statistical differences between algorithms. The symbol ‘+’ means that the algorithm has statistical difference with others, with better average. The symbol ‘−’ means that the algorithm has statistical difference with others, with worst average.
Full size table

As seen in [12], where the model was implemented with 2 objective functions, all algorithms failed to find the Pareto Front for the last operating hours. In this work, we incorporate a third objective function. The results show that The Pareto Fronts found by 4 out of the 5 algorithms reveal that the addition of the third objective function caused an improvement in the results. In particular, the Special-Crowd-Distance boosted MESH (E2V2D1C2) showed statistical robustness indicating that this approach benefits MESH for treating the OMRS problem, making it competitive compared to the most used algorithms from the pymoo library.

5 Conclusion

In this work, the performance of different evolutionary metaheuristics in solving a dispatch operation problem in cascade hydroelectric plants was analyzed. The work introduces a new version of the MESH algorithm, which added an alternative to its crowd distance operator, Special Crowding Distance. The mathematical modeling of the electrical problem was improved with the introduction of a third objective function. This function indicates the distance between the optimized water flow and the water flow in joint control mode. It guarantees a smoothing of the set of solutions found. Therefore, the work shows the first attempt to solve three conflicting objectives with the standard version and the Special-Crowd-Distance boosted MESH version. This version of MESH proved to be competitive compared to classic algorithms such as NSGA-II, NSGA-III, and MOEA/D, overlapping their results. This indicates that the Special-Crowd-Distance operation can benefit multi-objective algorithms based on Pareto Front. As future outcomes we indicate carry out a fine-tuning study of parameters for the new version of the MESH algorithm presented here, and also a decision-making study of which are the most suitable solutions for implementing electrical dispatch, taking into account the three objectives studied.