Applying the new multi-objective algorithms for the operation of a multi-reservoir system in hydropower plants (2024)

There are four significant storage dams (Gotvand Olia, Karun 1, Karun 3, and Karun 4) and a hydropower dam (Masjed-e-Soleyman) on the Karun River in the Karun basin, one of the richest basins in Iran: Due to 32 power plant units with a combined capacity of over 8500 MW (Table 1) for each of the dams, the system has a significant impact on supplying hydroelectric energy. With a nominal annual electricity output of 18,248.7 GW, it provides more than 90% of the nation's total hydroelectric energy (Dezab Consulting Engineering Co., 2019). The 67,257km² Karun watershed is located between 48° 10′ and 52° 30′ east longitude and 30° 20′ and 34° 05′ north latitude. Figure1 depicts the Gotvand Olia, Masjed-e-Soleyman, Karun 1, Karun 3, and Karun 4 dams, which are all situated on the main branch of the Karun River in succession, along with the location of the Karun basin.

Each hydropower plant comprises several independent power production units with a specific production capability (PPC). Each hydroelectric power facility also has a specific production efficiency to indicate its effectiveness. The quantity of electricity generated by each power plant is based on the efficiency coefficient over a particular time frame. The percentage of times a power plant runs at maximum efficiency throughout an operational period is denoted by this coefficient. The Karun basin’s electric reservoir system has proven successful over several months of operation, and each power plant can produce energy for at least 5–8h every day for at least 1month. Table 1 lists information about the electric power plants in the Karun River basin investigated for this study³⁵.

The specified system’s quantity and quality of exploitation are subject to stochastic influences, including the input flow to each reservoirs and multiple uses in the main electricity sector production of Iran, all of which multiply the importance of planning and justify any effort to optimize the objectives of the mentioned system, including hydropower generation.

Utilized data

In this study, 180-month time-series data, from 2000 to 2015, were used for the optimal operation of water resources in the Karun basin. The utilized data include net evaporation, river inflow, and total demands in the Gotvand Olia, Masjed-e-Soleyman, Karun 1, Karun 3, and Karun 4 dams (Table 2).

Problems with multi-objective optimization

While constrained by an MOP, multiple-objective functions must be minimized (or maximized) concurrently³⁶. Therefore, MOP can be explained as follows:

where \(f_{i} \left( {\vec{x}} \right)\) is the i-th singular objective function, R^d is the d-dimensional solution space, g_i and h_i are the i-th equal and unequal constraints, respectively.

In single-objective optimization, it is simple to compare various solutions using interface functions, and algorithms can only present one solution as the overall best solution. However, the interface operators lose their effectiveness in multi-objective optimizations because it is challenging to evaluate solutions with many objectives as the objectives contend with one another. It is essential to locate another interface operator for comparison with all objectives. In this instance, Pareto's dominance of the topic resulted in sensible solutions that weigh several objectives. The key definitions provided by Pareto for MOP minimization are listed below.

Definition 1

(Pareto dominance³⁷) vector \(\vec{x}\) dominates vector \(\vec{y}\) (denoted as \(\vec{x} \le \vec{y}\)) if:

Definition 2

(Pareto optimal) If and only if a strategy \(\vec{x}\) is considered Pareto optimal:

Definition 3

(Pareto optimal set) The definition of the Pareto optimum set (PS) is as follows:

Definition 4

(Pareto optimal front) This is how the Pareto optimum front (PF) is described:

When dealing with MOPs, finding a PS rather than a singular optimal solution is an obvious requirement. A typical technique for storing and recovering PS in the optimizer process creates an external archive with a particular number. Figure2 illustrates the relationship between the MOPs’ main ratios. In Fig.2, the continuous line represents PF, and x₁ is the PF's Pareto-optimal solution. x₁ is better than x₂, x₃, and x₄. Since x₃ has a higher value for f₂ and a smaller value for f₁ than x₄, x₃ and x₄ do not take precedence over each other.

Standard deviation (distance from the mean)

The standard deviation calculates the density of non-dominant solutions near the initial solution. MOPs are commonly used to maintain the distribution of solutions obtained in PF⁹. Based on the fitness values of each objective function, all of the solutions in the set of non-dominant solutions are arranged in ascending order. The average length of the rectangle sides for each solution, x_i, is determined by the neighboring solutions x_i−1 and x_i+1. These numbers stand for the solution xi's mean standard deviation. It is typical to assume that the distance from the set of solutions with the minimum and maximum values of the objective function is infinite. A greater range of solutions is possible at a greater distance from the density of solutions. Figure3 displays the separation from data congestion with two objectives.

Random archiving³⁸, grid archiving¹⁰, and ε-Pareto archiver³⁹ are some popular techniques for maintaining external archives. The random archiver randomly selects and removes the non-dominated solutions from the archive. The grid archiver partitions the target space into numerous hypercubes, and each hypercube saves the non-dominated solutions. Then, a preset number of the solutions are eliminated to distribute the remaining ones equally. A hypercube with two or more non-dominated solutions is archived using the ε-Pareto algorithm, but only the solution closest to the hypercube origin is kept. The value of ε determines how many hypercubes are broken up into the objective space.

Artificial hummingbird algorithm (AHA)

Zhao et al.¹³ developed the bio-based algorithm known as AHA to handle single-objective optimization problems. It can be categorized as a swarm algorithm. This program was created based on how hummingbirds naturally forage. Three foraging behaviors—guided, territorial, and migratory—are modeled by AHA during the optimization process. Axial, diagonal, and omnidirectional are the three types of flight mimicked by foraging behaviors. A visit table also imitates how hummingbirds use their memories to select their preferred food sources. The key components in AHA are briefly described below¹³.

Multi-objective artificial hummingbird algorithm (MOAHA)

MOAHA was first proposed by Zhao et al.⁴⁰ for MOPs. Three major components are integrated into MOAHA as follows to perform multi-objective optimization effectively:

1. (1)

   An external archive is introduced into MOAHA to save a fixed number of optimal non-dominated solutions during iterations.

2. (2)

   A dynamic elimination-based crowding distance (DECD) method is employed to efficiently maintain the archive during the entire optimization process, aiming to improve the solution diversity significantly.

3. (3)

   A solution update mechanism based on NDS is developed and used to refine non-dominated solutions in the solution update phase, remarkably contributing to solutions’ convergence.

Solution update strategy based on NDS

In AHA, in each iteration, after performing guided foraging or territorial foraging, the current solution is replaced with the candidate solution if the fitness value of the candidate solution is better than that of the current solution, or the current solution will remain unchanged. However, in MOAHA, dominance relation with many objectives makes solution updating complicated. The selection of solution update strategy will directly affect convergence towards the true PF.

In MOAHA, a solution update strategy based on NDS is proposed. In this strategy, all solutions are sorted by the NDS approach; solutions dominated by fewer solutions have a higher dominance hierarchy, and all solutions are assigned different ranks according to their non-dominated levels.

1. (1)

   When the front on which the candidate solution of the ith individual lies is better than the front on which the current solution of the ith individual lies, the current solution of the ith individual is replaced by its candidate solution. In Fig.4a, the candidate solution belongs to F₂ front, the current solution belongs to F₃ front, and the F₂ front is better than the F₃ front. In this case, the candidate solution will replace the current solution.

2. (2)

   When both the candidate and current solutions of the ith individual belong to the same front, there is a probability of 50% to choose between the two to update the current solution of the ith individual with the candidate solution. Figure4b depicts this case.

3. (3)

   When the front on which the candidate solution of the ith individual lies is worse than the front on which the current solution of the ith individual lies, the current solution will be retained for the next iteration. This case is described in Fig.4c. This solution update strategy is defined as follows:

   $$x_{i} \left( {t + 1} \right) = \left\{ {\begin{array}{*{20}l} {v_{i} \left( {t + 1} \right)} \hfill & {v_{i} \in F_{p} , x_{i} \in F_{q} \;and\;p < q} \hfill \\ {\left\{ {\begin{array}{*{20}l} {v_{i} \left( {t + 1} \right)} \hfill \\ {x_{i} \left( t \right)} \hfill \\ \end{array} } \right.} \hfill & {\begin{array}{*{20}l} {v_{i} \in F_{p} , x_{i} \in F_{q} ,\quad p = q\;and\;rand < 0.5} \hfill \\ {v_{i} \in F_{p} , x_{i} \in F_{q} ,\quad p = q\;and\;rand \ge 0.5} \hfill \\ \end{array} } \hfill \\ {x_{i} \left( t \right) } \hfill & {v_{i} \in F_{p} , x_{i} \in F_{q} ,\quad and\;p > q} \hfill \\ \end{array} } \right.$$

   (8)

where F_p is the pth front in NDS.

Dynamic elimination-based crowding distance (DECD)

The crowding distance is an effective parameter-free method to enhance the diversity of solutions⁹, so it is frequently employed in many MOAs to maintain a fixed-size external archive by removing excessive solutions with smaller crowding distance. However, whenever a solution is removed, the sorting of the crowding distance of the remaining solutions will be changed. Then, the next solution to be removed in terms of the initial crowding sorting may not be the one with current crowding sorting, which will reduce solution diversity in the external archive. Figure5 displays the sorting change of the crowding distance before and after a solution is removed. Evidently, sorting the subsequent crowding distance is changed when the solution with the smallest crowding distance (x₄) is removed. The next solution to be removed will be x₅ according to the original sorting of the crowding distance in Fig.5a.

Nevertheless, x₆ should be removed according to the updated sorting of crowding distance in Fig.5b. The archiving procedure based on crowding distance is ineffective in maintaining the sustainable diversity of PS. Thus, the external archive with DECD is proposed. When the solution with the current minimal crowding distance in PS is removed, only the crowding distances of the immediately adjacent solutions to remove the solution must be updated, and the crowding distances of the rest of the solutions do not change.

Before the ith solution with the current smallest crowding distance is removed, the crowding distances of the last and next solutions adjacent to the ith solution are given as follows, respectively:

After the ith solution with the current smallest crowding distance is removed, the crowding distances of the last and next solutions adjacent to the ith solution are changed as follows, respectively:

where \(d_{k} \left( {x_{i - 1} } \right)\) is the side length of the cuboid of the (i−1)th solution, and \({\text{D}}\left( {x_{i - 1} } \right)\) is the crowding distance of the (i−1)th solution. Figure6 illustrates the change of the crowding distance of the solutions adjacent to the removed one using DECD for two objectives. When n solutions are removed for establishing an external archive using the DECD method, only crowding distances of 2n solutions need to be recalculated. Liu and Chen⁴¹ proposed a crowding distance elimination (CDE) method for improving solution diversity. Luo et al.⁴² devised a dynamic crowding distance-based diversity maintenance strategy (DCDDMS) to enhance solution diversity horizontally. DECD, CDE, and DCD-DMS are similar in that the sorting of the residual solutions is reconsidered whenever a solution with the minimum crowding distance so far is removed. The difference among them is that whenever the solution with the current minimum crowding distance is removed, the crowding distances of residual solutions in the other methods are recalculated. In contrast, only the crowding distances of the immediately adjacent solutions to the removed solution in the DECD method are recalculated. The DECD method reduces the computational complexity of computing the crowding distance of solutions from O(TMN2) to O(TMN).

Figure7 shows external archiving procedures based on crowding distance and DECD for two objectives. There are 100 non-dominated solutions in Fig.7a, in which 50 must be removed. Figure7b, c display the external archives based on crowding distance and DECD, respectively. Figure8 shows external archiving procedures based on crowding distance and DECD for three objectives. There are 100 non-dominated solutions in Fig.8a, in which 50 must be removed. Figure8b, c depict external archives based on crowding distance and DECD, respectively. Evidently, the solutions in the external archiving based on the DECD method maintain better diversity and distribution than those in the external archiving based on the crowding distance.

Guided foraging

In AHA, a hummingbird tends to move to the food source with the highest visit level when performing guided foraging. When there are multiple food sources with the highest visit level, the one with the best fitness is selected as the target food source for foraging. In MOAHA, when there are multiple food sources with the highest visit level, the one in the better front based on NDS is selected as the target food source. If there are still multiple solutions in the same better front, one is randomly selected as the target food source.

After a hummingbird visits the target food source using Eq.(8), the visit level of the target food source for this hummingbird in the visit table is set to 0, and then the visit levels of the other food sources are increased by 1. After performing guided foraging, a new food source (candidate solution) is obtained, and the solution update strategy based on NDS is performed.

If this new food source meets the requirement of update condition in Eq.(9), this hummingbird will abandon the previous food source and station in this new food source (candidate solution). Otherwise, it still rests on the previous food source. This solution update results in the visit level update of the food source for all other hummingbirds. If the solution update succeeds for a hummingbird, the visit level of the updated food source should be set to the highest level of the other food sources plus 1 for every other hummingbird. If the solution update fails for a hummingbird, the visit level does not change.

The following example (a two-objective minimization problem) shows how each hummingbird's target food source is selected and how the visit table is modified when employing the solution update strategy.

In Fig.9a, assume that four hummingbirds are placed on different food sources, and all visit levels in the visit table are initialized to 0. The first hummingbird finds the food sources of other hummingbirds x₂, x_3, and x₄ that have the same highest level; since x₃ dominates both x₂ and x₄, the food source of the hummingbird x₃ is treated as the target source of the first hummingbird. After the first hummingbird performs guided foraging via Eq.(8), the visit level of the food source of the hummingbird x₃ is set to 0, and the visit levels of the food sources of hummingbirds x₂ and x₄ are increased by 1. Since the current solution × 1 dominates the candidate solution v₁, x₁ is not updated by v₁ (if a solution dominates another solution, the front of this solution is better than the front of the other solution).

In Fig.9b, the second hummingbird finds the food sources of other hummingbirds x₁, x₃ and x₄ with the same highest level, since x₁ dominates both x₃ and x₄, the food source of hummingbird x₁ is treated as the target source of the second hummingbird. After the second hummingbird performs guided foraging via Eq.(8), the visit level of the food source of the hummingbird x₁ is set to 0, and visit levels of the food sources of the hummingbirds x₂ and x₄ are increased by 1. Since the candidate solution v₂ dominates the current solution x₂, x₂ is updated by v_2, and the visit level of the source of the hummingbird x₂ for each of the other hummingbirds is modified to the highest visit level increased by 1 in every corresponding row.

In Fig.9c, the third hummingbird finds the food sources of the hummingbird x₂ that has the highest level, so the food source of the hummingbird x₂ is treated as the target source of the third hummingbird. After the third hummingbird performs guided foraging via Eq.(8), the visit level of the food source of the hummingbird x₂ is set to 0, and visit levels of food sources of the hummingbirds x₁ and x₄ are increased by 1. Since the current solution x₃ dominates the candidate solution v₃, x₃ is not updated by v₃.

In Fig.9d, the fourth hummingbird finds the food source of the hummingbird × 2 that has the highest visit level, and the food source of the hummingbird × 2 is treated as the target source of the fourth hummingbird. After the fourth hummingbird performs guided foraging via Eq.(8), the visit level of the food source of the hummingbird × 2 is set to 0, and visit levels of the food sources of the hummingbirds × 1 and × 3 are increased by 1. Since candidate solution v4 dominates the current solution × 4, × 4 is updated by v4, and the visit level of the source of the hummingbird × 4 for each of the other hummingbirds is modified to the highest visit level increased by 1 in every corresponding row. When one iteration is achieved for the four hummingbirds, their updated solutions and visit table are shown in Fig.10.

Territorial foraging

When a hummingbird performs a local search in its territory, a new food source is produced as a candidate solution. The mathematic equation simulating the territorial foraging of hummingbirds is given below to improve the solution diversity in MOAHA:

where xa(t) is a solution randomly chosen from the external archive.

When a hummingbird performs territorial foraging, the visit levels of the food sources of other hummingbirds are increased by 1. After completing the territorial foraging, a new food source (candidate solution) is found, and the solution update is performed according to Eq.(13). If the solution update succeeds, the visit level of the updated food source is set to the highest level of the other food sources plus 1; if the solution update fails, the visit level does not change. The pseudocode of the territorial foraging strategy of MOAHA is described in Fig.9, showing the update of the visit table when the first hummingbird performs territorial foraging.

For brevity, Fig.11 only shows the solution update procedure of the first hummingbird in territorial foraging. Based on Fig.11a, when the first hummingbird performs territorial foraging in its local region, the visit levels of other hummingbirds × 2, × 3, and × 4 are increased by 1. After performing territorial foraging, since candidate solution v1 dominates the previous solution × 1, the previous solution × 1 is replaced by the candidate solution v1, and the source of the hummingbird × 1 for each of the other hummingbirds is modified to the highest visit level increased by 1 in every corresponding row. Figure11b shows the updated solution and visit table after the first hummingbird performs the territorial foraging.

Note that the solution update procedure is the same for both guided foraging and territorial foraging; the difference between them in the updating visit table is that the visit level of the visited food source is initialized to 0 in guided foraging, while this step is not implemented in territorial foraging as a hummingbird visits its own territory instead of the existing food source of other hummingbirds.

Migration foraging

For AHA, when performing migration foraging, the hummingbird located at the food source with the worst nectar-refilling rate tends to migrate to a new food source produced randomly in the entire search space. The food source with the worst nectar-refilling rate is the solution with the worst fitness in all the solutions. For MOAHA, the food sources with the worst nectar-refilling rate are defined as the solutions in the worst front based on NDS.

The migration foraging of hummingbirds in MOAHA is expressed as follows:

where Fend is the worst front, r is the rand number in [0, 1], and Up and Low are the upper and lower boundaries. When migration foraging is completed, the visit table is updated, and the update procedure is similar to those of the other two foraging strategies.

Multi-objective artificial hummingbird algorithm (MOAHA)

MOAHA begins with a random group of solutions by creating an external archive with a fixed number for a particular multi-objective optimization task and initializing a random population of hummingbirds. After initializing a visit table, all non-dominated solutions from the original population are recorded in the archive. Each time, MOAHA has a 50% chance of choosing between territorial foraging and guided foraging. When a hummingbird participates in guided foraging, it usually shifts its position concerning the target food supply determined by the visit table and dominance relationship. When a hummingbird is territorially foraging, it shifts its position toward the community in its area. The solution update based on NDS is performed in line with Eq.(8) following the completion of one forage, which modifies the visit table.

Migration foraging is used every 2n iterations, and the visit table and worst-front solutions are both randomly initialized in the search field.

At the end of each iteration, non-dominated solutions in the new population are added to the archive. The external archiving mechanism based on DECD is triggered if the archive size exceeds the predetermined limit. These processes are repeated until the maximum number of iterations is reached. The PF is eventually returned along with the ideal non-dominated solutions from the archive. The MOAHA diagram is explained in Fig.12.

This is a theoretical investigation of combining different techniques.

1. (1)

   MOAHA has an internal external archive with a preset size that can store and recover non-dominated optimal solutions. Archives may be a great resource for all MOAs. Solutions are managed in archives using the DECD technique. The literature has established that the crowing distance method is suitable for maintaining solution variety. Contrarily, the one-step elimination method of crowding distance is useful for maintaining archives. The phase-out technique can greatly improve the distribution of solutions. As a result, the DECD approach can considerably improve the consistency and diversity of the solutions.

2. (2)

   The NDS approach can divide all the solutions into multiple fronts, allowing for easy comparison of both the existing and potential solutions based on front ranks⁹. Better non-dominated solutions can then be archived and handed on to the following iteration to aid in the search operations of the other users. Therefore, a solution updating method with NDS can improve the algorithm's ability to converge to optimal solutions.

3. (3)

   Since MOAHA does a global search using the same methods as AHA, convergence may result. The performance of searches in terms of exploration, exploitation, and convergence has been theoretically assessed for these procedures, three foraging methods, and a visit table⁴⁰. Because MOAHA receives all the advantages of AHA, searchers can perform both exploitation and exploration in the same manner. The main difference between the two is that while AHA simply stores one optimal solution and later fine-tunes it, MOAHA looks through various non-dominated solutions and saves them to the archives.

Other multi-objective algorithms

Multi-objective moth swarm algorithm (MOMSA)

Single-objective MSA, proposed by Mohamed et al.⁴³, is inspired by the behavior of moths in nature. Moths try to hide from predators during the day while looking for food resources at night with a celestial navigation technique. They fly in a straight line over a long distance by steering their locomotion in a steady angle relative to moonlight as the celestial far-distant point of light. In the MSA, the possible solution is represented by the position of the light source, and the quality of this solution is considered the luminescence intensity of the light source. Three groups of moths (pathfinders, prospectors, and onlookers) are considered in the MSA. Pathfinders can find the best position over the optimization space with the First-In, Last-Out principle to guide the movement of the main swarm. Prospectors tend to wander into a random spiral path near the light sources marked by the pathfinders. Onlookers drift directly toward the best global solution (moonlight) achieved by the prospectors⁴³.

The dominant features of the algorithm (such as moonlight and pathfinder moth) must be properly defined to convert MSA into an efficient multi-objective optimization algorithm. In normal optimization problems for MSA, only one objective function should be minimized or maximized. In this condition, three types of moths, including pathfinder, prospector, and onlooker, follow the process of searching and reaching the final objective in three stages of cognition, transverse orientation, and constellation orientation. The prospector moths randomly move in a maze toward the best pathfinder moths (responsible for guiding the prospectors). Following this process, the onlooker moths move directly toward the best global optimal solution (moonlight) detected by the prospectors. The extraordinary performance of the MSA algorithm in solving single-objective problems, proven in various studies, has turned it into one of the most powerful metaheuristic algorithms^44,45. A multi-objective moth swarm algorithm (MOMSA) with the capability of solving complex and large-scale problems was proposed by Sharifi et al.⁴⁴ to solve multi-objective constrained problems. Since more than one function must be maximized in MOPs, the definition for selecting the type of moths and the best value (moonlight) is changed in MSA to the multi-objective space. The crowding-distance mechanism selects the most efficient (best) solutions in the population as the pathfinder moths and moonlight (global optimum). Sharifi et al.⁴⁴ presented a full description of the method and key components of the MOMSA algorithm.

Multi-objective material generation algorithm (MOMGA)

The MGA algorithm, a metaheuristic algorithm recently developed by Talatahari et al.⁴⁶, considers the general aspects of material generation in nature, along with the basic and advanced principles of chemistry, including chemical compounds and chemical reactions⁴⁶.

The multi-objective version of MGA for solving multi-objective optimization problems was presented by Nouhi et al.⁴⁷. MOMGA has three new mechanisms for solving multi-objective optimization.

The first mechanism introduced into MGA is the archive, which serves as a storage facility for storing or restoring the derived Pareto optimal solutions. The archive has a single controller that manages which solutions are added to the archive and when the archive is full. There is a limit to the number of solutions that can be stored in the archive. The residents of the archive are compared to the non-dominated solutions created. Three major scenarios are possible:

1. I.

   If at least one member in the archive dominates the new solution, it is not allowed to enter the archive.

2. II.

   The new solution may be added to the archive if it dominates at least one solution in the archive by omitting the one already in the archive.

3. III.

   If the new and archive solutions do not dominate each other, the new solution is added to the archive.

The grid mechanism included in MGA is the second effective technique for enhancing non-dominated solutions in the archive. If the archive is full, the grid technique will be utilized to reorganize the object space's segmentation and find the most populated area to eliminate a solution. The additional member should be included in the least crowded segment to boost the variety of the final approximated Pareto optimal front. As the number of possible solutions in the hypercube increases, the possibility of deleting a solution rises. If the archive is full, the most crowded areas are chosen first, and a solution from one of them is randomly deleted to make way for the new solution. When a solution is placed outside the hypercubes, a special case arises. All segments in this scenario have been expanded to fit the most recent solutions. As a result, the segments of alternative solutions can also be changed.

Due to Pareto optimality, solutions in a multi-objective search space cannot be compared; hence, the leader selection mechanism is the last machine in MGA. Consequently, MGA includes a leader selection method to address this problem.

The search leaders guide the other search candidates to possible areas of the search space to attain a solution close to the global optimum. As previously stated, the archive contains only the best non-dominated solutions. The leader selection mechanism chooses the least crowded portions of the search space and presents the best as non-dominated answers.

A full description of the MOMGA algorithm is provided by Nouhi et al.⁴⁷.

Performance criteria for MOPs

GD⁴⁸, S⁴⁹, Δ⁵⁰, and MS³⁶ are the four performance parameters used in this study to assess the performance of the exploitation algorithms.

where d is the Euclidean distance between the i-th member in PF_g and the closest member in PF_optimal, where NPF is the total number of members discovered in the Pareto front. The best GD criterion that can be obtained is zero, which indicates that PF_g is precisely on the optimal line or PF_optimal.

In Eq.(16), \(d_{i} = min_{j} (\left| {f_{1}^{i} \left( x \right) - f_{1}^{j} \left( x \right)} \right| + \left| {f_{2}^{i} \left( x \right) - f_{2}^{j} \left( x \right)} \right|\), i,j = 1,2,…,NPF and \(\overline{d}\) is the average of all d_i. The lowest value of S in PF_g leads to the best uniform distribution. As a result, if all non-dominated solutions are distributed equally throughout PF_g, the value of S is zero.

where in PF_optimal and PF_g, df and dl represent the separation between extremum solutions (beginning and finishing points). The distance di separate any point in PFg and the closest point in PFotpimal.

This metric represents a graph's point-to-point distance. The value of Δ is always greater than zero, and a lower value indicates a better distribution and spread of solutions.

When Δ is zero, it is a perfect condition showing that \(d_{i} = \overline{d}\) for all non-dominant points.

where \(f_{i}^{max}\) and \(f_{i}^{min}\) are, respectively, the maximum and minimum values of the ith objective in PF_g, and \(F_{i}^{max}\) and \(F_{i}^{min}\) are the maximum and minimum of the ith target in PF_optimal, respectively. A higher MS score indicates a wider range of solutions.

Generally, the performance criteria of multi-objective problems are used to evaluate the performance of multi-objective algorithms. The closer the GD, S, and Δ criteria to zero and the higher the MS criterion, the greater the superiority and performance of the multi-objective algorithm.

Long-term multi reservoir-multi objective simulation–optimization model

The optimal monthly release rates from the dams' reservoirs serve as the decision factors in the multi-reservoir, multi-objective optimization model. In succession, these release values in the system of five Karun dam reservoirs include releases intended to satisfy downstream demand, control flooding, and provide electric energy. The MOAHA algorithm has 1800 decision variables in the system because the planning horizon for this study is 180 months (from 2000 to 2015). Reservoir downstream needs and releases from energy outputs are the decision variables, while storage volume and input to reservoirs are the state variables for each time. The volume of the river, the height of evaporation, the height of precipitation, and the volume of requirements every month are the input data for the model.

For deterministic optimization of the reservoir system in reservoir problems, the usual objective function is as follows:

In this context, t stands for the intended period index, T for the number of exploitation periods, Z for the target that should be minimized or maximized, \(Re_{t}\) is the release, and \(S_{t}\) is the storage in period t. As a fundamental relationship in the state-space equations, the conservation of mass in the system is as follows:

Thus, \(Q_{t}\) is the inflow during period t, \(Sp_{t}\) is the overflow during period t, and \(Loss_{t}\) is the loss during period t. The upper and lower storage volume limitations should be set to meet the needs of renewable flood control volume objectives and ensure a minimum level of balance for reservoir dead volume and power plant operation.

In these equations, \(S_{Min}\) and \(S_{Max}\) are the minimum and maximum storage volumes,\(Re_{Min}\) and \(Re_{Max}\) are the minimum and maximum release volumes, respectively.

Evaporation is used to evaluate reservoir losses while accounting for the non-linear relationship between reservoir surface and volume based on Eqs.(23) to (25):

In these equations, a, b, and c are constant coefficients for converting reservoir volume to the corresponding level of the reservoir, respectively, \(Ev_{t}\) the average drop during period t (evaporation minus precipitation) in millimeters, \(\overline{{A_{t} }}\) the average level of reservoir during period t in square kilometers, and \(A_{t}\) and \(A_{t + 1}\) the amount of level in reservoir at beginning and end of period t.

The objectives of providing downstream needs, flood control, and electric energy production are considered for the optimal use of the reservoir system. Three outlets—an electrical outlet, a downstream demand supply outlet, and a flood control outlet—are considered for the reservoir system. In various exploitation modes, water is released first from the electrical outlet, then from the supply outlet, and finally, from the flood control outlet.

In this way, the other two outlets are not used until the electricity outlet is used to its maximum capacity, and after using the demand–supply outlet to its maximum capacity, the flood control outlet is utilized. \(Re_{Max}^{Power}\) is the maximum output capacity for power production, \(Re_{Max}^{De}\) is the maximum output capacity for demand supply, and \(Re_{Max}^{FC}\) is the maximum output capacity for flood control, all measured in million cubic meters. Note that in this instance, \(Re_{t}\) is between zero and the total of the maximum output of reservoir, and \(Re_{t}^{Power}\), \(Re_{t}^{De}\) and \(Re_{t}^{FC}\) respectively denote the quantity of release from electric output, demand supply, and flood control measured in a million cubic meters for each period t.

Operating the reservoir system for downstream demands supply

In this case, the goal is to minimize the sum of squared difference from the requirement in period t, displayed as the following equation:

Here, \(f_{1}\) is the sum of the squared difference of monthly release (\(Re_{t}^{Power} + Re_{t}^{De}\)) from requirement in period t, \(De_{Max}\) is the maximum monthly requirement of the reservoir, and \(De_{t}\) is the amount of requirement.

Operating the reservoir system for flood control

Here, all the equations are like the exploitation for satisfying downstream needs, and only the objective function is as follows:

Here, \(S_{t}^{target}\) is the desired flood control volume in t period per million cubic meters. In the introduced objective function, the objective is to keep a constant volume around \(S_{t}^{target}\) in all the exploitation periods. In fact, if the reservoir volume is more than the desired volume, the objective is flood control; if it is less than the desired volume, the objective of supplying downstream needs will not be met and will cause deviation from \(S_{t}^{target}\) in each period.

Operating the reservoir system for hydropower energy generation

One must consider the problem's complexity in terms of constraints and nonlinear conditions to make the best use of electricity. For this purpose, the final objective function, including minimizing the lack of production power compared to the installed capacity of power plants, was used according to Eq.(28).

Here, \(P_{i,t}\) is the production energy of reservoir i in period t per MW, and \(PPC_{i} { }\) is the installation capacity of reservoir i per MW. According to Eqs.(29) to (40), additional constraints and criteria were applied to optimize hydropower production:

where \({\text{e}}_{{{\text{i}},{\text{t}}}}\) is the reservoir power plant’s efficiency, assumed to be constant over all periods, and g is the acceleration of the earth’s gravity, which is equal to 9.81 m per square second. \({\text{PF}}_{{\text{i}}}\) is the reservoir power plant’s operating coefficient, and \({\text{Mul}}_{{\text{t}}}\) is the ratio of million cubic meters to cubic meters per second over the course of a time interval of t. \({\text{H}}_{{{\text{i}},{\text{t}}}}\) represents the average water level of reservoir i for t meters, \({\text{H}}_{{{\text{i}},{\text{t}} + 1}}\) represents the water level of reservoir i at the start of period t (meters), \(TW_{i,t}\) represents the water level of reservoir i at the end of period t (meters), \(a_{0i}\), \(a_{1i}\), \(a_{2i}\), \(a_{3i}\) and \(a_{4i}\) indicate the constant coefficients used to convert reservoir volume to the corresponding height in reservoir i, constant coefficients of conversion \({\text{b}}_{{0{\text{i}}}}\), \({\text{b}}_{{1{\text{i}}}}\), \({\text{b}}_{{2{\text{i}}}}\), \({\text{b}}_{{3{\text{i}}}}\) and \({\text{b}}_{{4{\text{i}}}}\) output water from the power plant to the elevation of the reservoir i, \({\text{RPS}}_{{{\text{i}},{\text{t}}}}\) dentoes the volume of water that overflows from reservoir i's hydropower outlet during the t period, and \({\text{RP}}_{{{\text{i}},{\text{t}}}}\) shows the volume of water released from reservoir i's hydroelectric outlet during period t to generate electricity. Reservoir i's natural river influx during period t is represented by \({\text{Q}}_{{{\text{i}},{\text{t}}}}\), its overflow is denoted by \({\text{Sp}}_{{{\text{i}},{\text{t}}}}\), and its overall volume loss over time is represented by \({\text{Loss}}_{{{\text{i}},{\text{t}}}}\). \({\text{A}}_{{{\text{i}},{\text{t}}}}\) and \({\text{A}}_{{{\text{i}},{\text{t}} + 1}}\), which shows the reservoir i's water level at the start and end of t period, respectively, represents the reservoir i's water level in square kilometers. \({\text{Ev}}_{{{\text{i}},{\text{t}}}}\) represents the amount of reservoir i's water lost in millimeters in period t (evaporation minus precipitation). \({\text{S}}_{{{\text{Min}}\;{\text{i}}}}\) is the minimum storage volume of the reservoir i, \({\text{S}}_{{{\text{Max}}\;{\text{i}}}}\) is the maximum storage volume of the reservoir i and \({\text{Re}}_{{{\text{Min}}\;{\text{i}}}}\) and \({\text{Re}}_{{{\text{Max}}\;{\text{i}}}}\) are the minimum and maximum release volumes of the reservoir i, respectively. The constant coefficients of conversion for the reservoir volume to the corresponding level of the reservoir are \({\text{c}}_{{0{\text{i}}}}\), \({\text{c}}_{{1{\text{i}}}}\), \({\text{c}}_{{2{\text{i}}}}\), \({\text{c}}_{{3{\text{i}}}}\) and \({\text{c}}_{{4{\text{i}}}}\). The reservoir i's lowest and maximum release volumes are denoted by \({\text{Re}}_{{{\text{Min}}\;{\text{i}}}}\) and \({\text{Re}}_{{{\text{Max}}\;{\text{i}}}}\), respectively. \({\text{Re}}_{{{\text{i}},{\text{t}}}}^{{{\text{Power}}}}\) is a function of \({\text{Re}}_{{{\text{i}} - 1,{\text{t}}}}^{{{\text{Power}}}}\). The power and amount of hydroelectric energy the reservoir discharged over time t. Power is the amount of water released from a reservoir's hydroelectric output over a given time span t.

Performance metrics for reservoirs

The final and most crucial phase in employing optimization and simulation models to exploit reservoirs is the evaluation of operation policies.

Reliability

Reliability can be defined in two ways: time and volume. Time reliability means the percentage of periods when the system fully supplies the existing needs and does not fail. The following equation is used to determine this parameter’s value.

where NDe_f stands for the overall number of failures experienced during the operation period, Dei for the demand value experienced during the i period, Re_i for the output value experienced during the i period, and Rel for the system's reliability experienced during the i period. The higher the value of this parameter, the higher the time reliability of the system⁵¹.

Vulnerability

This indicator measures the intensity of system failures and is derived from the equation below⁵¹:

where Val is the system’s vulnerability, De_i is the amount of demand in period i, Re_i is the output value in period I, and t is the total number of exploitation periods.

Resiliency

This index demonstrates the system’s capacity to alter the state of affairs. Reversibility is the likelihood that a system will resume its desired condition after failure⁵¹. The equation below can be used to calculate this parameter’s value.

where Def_i is the deficiency in the i-th period, and \(\begin{array}{*{20}c} T \\ N \\ {i = 1} \\ \end{array} ()\) is the number of times the condition in parentheses happened.

Sustainability index

This index aggregates system performance requirements into a general index based on the following reservoir performance indicators, and aims to enable comparison and decision-making between various scenarios⁵².

Reliability, resiliency, vulnerability, and sustainability criteria are used to compare different algorithms in reservoir operation problems. The closer the reliability, resiliency, and sustainability indicators to 100% and the closer the vulnerability index to 0%, the greater the success in the investigated problem.

Multi-objective criteria problems

Problem 1 of the Schaffer benchmark

A benchmark problem was provided by Schaffer⁵³ to assess the effectiveness of multi-objective optimization techniques. This univariate convex problem has two objectives that must be minimized (Eq.45). This test problem has served as a reliable benchmark. Almost every multi-objective EA has been tested on this problem.

This criterion seeks to minimize f1 and f2 simultaneously. The Pareto optimum sites are located in the region where x ∈ [0, 2], as shown in Fig. 13. There is a trade-off between two functions that grow and decrease within this interval because f₁ and f₂ increase outside of it.

Problem 2 of the MMF1 benchmark

A benchmark concern for CEC 2020 is MMF1. The most recent and official multimodal-multiobjective benchmark set is the CEC 2020. The CEC 2020 benchmark set includes four distinct kinds of benchmark problems: convex linear, convex nonlinear, concave, and concave linear. The shapes of the objective functions and the spaces of the choice variables distinguish these benchmark problems.

The Pareto front (PF) and Pareto optimum set (PS) equations and forms for MMF1 are as follows:

where the search space is x₁ ∈ [− 1, 1] and x₂ ∈ [1, 3].

Pareto global sets (PSs) are:

where \(1 \le x_{1} \le 3\) , and the global Pareto Front (PF) is:

where \(0 \le f_{1} \le 1\). Its actual PS and PF are shown in Fig.14.

Liang et al.⁵⁴ provide more details on the additional characteristics of the CEC benchmark set problems. The MOAHA multi-objective algorithm's performance in this work was tested and verified using the two benchmark challenges described above.