Deriving Operating Rules of Hydropower Reservoirs Using Multi-Strategy Ensemble Henry Gas Solubility Optimization-Driven Support Vector Machine

Abstract

Hydropower is an important clean renewable energy that plays a key role in coping with issues such as global energy security, environmental protection, and climate change. In order to improve the optimal operation ability of hydropower reservoirs in the context of forecast runoff with limited accuracy and prediction period, there has been a growing interest in deriving operating rules of hydropower reservoirs. Reasonable operation decision is very important for safe operation of reservoirs and efficient utilization of water resources. Therefore, a novel method of operation rules derivation is proposed in this study. Optimal operation model of hydropower reservoir is established and support vector machine (SVM) is used to derive operation rules based on the optimal operation results. In order to improve the performance of SVM, the Henry gas solubility optimization (HGSO) is used to optimize its hyperparameters for the first time. Meanwhile, multiple strategies are applied to overcome the drawbacks of HGSO. The multi-verse optimizer (MVO) is used to enhance the exploration capability of basic HGSO. Quadratic interpolation (QI) is used to improve the exploitation ability of HGSO. In this study, the Xiluodu and Xiangjiaba hydropower reservoirs in the upper Yangtze River of China were selected as a case study. First, the improved HGSO called MVQIHGSO was tested on 23 classical benchmark functions. Then, it was employed to optimize hyperparameters of SVM model for deriving operation rules. The results and statistical studies indicate that the improved HGSO outperforms the comparison algorithms in exploration and exploitation. The obtained results imply that the novel method named MVQIHGSO-SVM can provide a new practical tool to deriving operation rules for hydropower reservoirs, which is conducive to the safe and efficient utilization of water resources.

1. Introduction

Climate change has undoubtedly emerged as the great challenging global problem facing humanity in this era, which significantly threatens the human life, industrial production, and ecological environment [1,2,3,4]. The main reason of increased climate change is the large amount and rapid growth of non-renewable energy consumption over the past several decades, mainly fossil fuels [5,6]. In response to climate change, global warming will be limited to “well below 2 °C” before the level of industrialization, and efforts will be made to control it at 1.5 °C [7]. Meanwhile, clean renewable resources including hydropower, solar energy, wind energy, geothermal energy, biomass energy, and tidal power have attracted widespread attention from policy-makers and academics in all developed and developing countries of the world to deal with the problem of ecological deficit [8,9]. Particularly, hydropower is the most widely used renewable energy in human society due to its advantages of low operating cost, stability, and flexible on demand power supply as well as other comprehensive benefits such as flood control, irrigation, navigation, and ecological environment protection [10,11]. However, the hydropower utilization is inevitably limited by natural conditions such as hydrology, climate, and geomorphology [12]. Actually, to solve this problem, reservoirs’ optimal operation is able to change the temporal and spatial distribution of water resources, which is an important mean and tool to realize the rational allocation and efficient utilization of water resources [13].

Generally, the mathematical programming techniques and meta-heuristic optimization methods are the two approaches used for optimal operation of reservoirs based on the inflow runoff [14,15,16,17]. Specifically, the mathematical programming can be divided into three types: linear programming (LP), non-linear programming (NLP), and dynamic programming (DP) as well as its many variants such as incremental dynamic programming (IDP), successive approximation DP (DPSA), discrete differential DP (DDDP), and progressive optimality algorithm (POA) [18]. LP is the most widely used optimization technology in the field of water resources planning and management due to its simplicity and adaptability; however, it cannot deal with the non-linear constraints of the reservoir optimal operation [19]. To solve this problem, NLP is applied to optimize the reservoir storage and release process; however, the application of NLP to large scale reservoirs optimal operation is limited by the high computational complexity [19,20]. The aim of DP is to overcome the inherent problems of reservoir operation including non-linear and stochastic features, which has been the most popular optimization method so far. Nevertheless, it would suffer from “curve of dimensionality” with the reservoir scales increase [19,21]. Fortunately, many variants of DP mainly composed of IDP, DPSA, DDDP, and POA were developed to overcome the dimensionality problem in the past decades [22]. On the other hand, the emergence of meta-heuristic optimization methods provides a novel and promising way for solving reservoir optimal operation model, which only requires objective function information without derivatives or other auxiliary knowledge [23]. Although a large number of research studies have focused on the application of these method to reservoir operation, they are rarely applied in the actual reservoir engineering operation due to the randomness of results and low convergence speed. In addition, the above techniques may lead to ineffective and inefficient operation under the changing environment conditions. For example, the operator needs to carry out repeated trial calculations to formulate a better reservoir operation scheme because of the uncertainties in reservoir operations, leading to great difficulties to the work of decision makers.

In order to formulate reasonable operation scheme with high efficiency, a simple, scalable, and adaptable operation rule is important and necessary for operators [24]. Operation diagram and operation function are the two types of operation rules [25]. The former is a control curve diagram to guide the operation of the reservoir by dividing the storage capacity of the reservoir into different operation area through some index line formed by water storage and water supply. The operation diagram is characterized by simplicity, intuition, and easy operation; however, the application of operation diagram is often affected by the subjectivity of operators. Therefore, there is still ample room to improve the economic benefits of the reservoirs. Fortunately, operation function is another form of operation rule with clearly defined input and output to make operational decisions considering uncertain conditions [26]. For example, the end water level as decision variable during a period can be determined by certain input, such as initial water level, inflow runoff, and other hydrometeorological information. This method has been regarded as the most promising future lines of research in this field [27,28,29]. The operation rule derivation methods based on operation function can be divided into two types: regression analysis and machine learning technique [30]. The common point of these two methods is that the optimal operation process of reservoir groups should be obtained according to the implicit stochastic optimization method [31]. The significant difference between the two method is that the expression of regression analysis is based on speculation, inspection and correction, in which the process is complex and greatly affected by human factors; however, the latter can solve this problem well and obtains the reservoir group operation function with good adaptability.

Support vector machine (SVM) is a newly developed machine learning method based on statistical learning [32,33]. It has the advantages of global optimization and good generalization performance, and the diversity of kernel functions also solves different types of data simulation and classification problem in processing small sample data. As we all know, hyperparameter setting has a direct impact on the performance of SVM. At present, grid search method and meta-heuristic intelligence algorithms are two commonly used methods to optimize hyperparameter. In essence, the former is an enumeration procedure used to determine the combinations of available hyperparameters, which is very time-consuming [34]. The latter assumes that a good combination of hyperparameters has an aggregation effect, which is able to find the most suitable hyperparameter combination in a shorter time [35]. Therefore, the hybrid model combined meta-heuristic optimization algorithm and machine learning method is applied for the operation rule derivation of reservoir groups in this study. In order to improve the derivation performance, a new meta-heuristic optimization algorithm namely Henry gas solubility optimization, short for HGSO, is introduced to optimize the hyperparameters of SVM for the first time. Unfortunately, it is found that the basic HGSO algorithm has some disadvantages such as local optima stagnation and low convergence speed [36]. To handle these drawbacks of the basic HGSO algorithm, multiple strategies with multi-verse optimizer (MVO) and quadratic interpolation method (QI) are firstly applied to enhance the exploration and exploitation ability of basic HGSO, called MVQIHGSO algorithm in the present paper. As a result, the hybrid model combined MVQIHGSO algorithm and SVM method, called MVQIHGSO-SVM is proposed to derivate operation rule of hydropower reservoirs to maximize the hydropower benefits.

A novel operation rule derivation method is proposed based on the improved HGSO algorithm and SVM. The reminder of this paper is organized as follows. In Section 2, the process of operation rule derivation using the proposed MVQIHGSO-SVM method was introduced, including the optimal operation model of hydropower reservoirs, the basic HGSO algorithm and its improvement strategies as well as the basic principle of SVM. In Section 3, numerical experiments were carried out to validate the performance of the MVQIHGSO algorithm. In Section 4, the Xiluodu and Xiangjiaba hydropower reservoirs in the lower reaches of the Jinsha river of the upper Yangtze river in China were selected as a case study. In Section 5, the conclusions were drawn.

2. Materials and Methods

2.1. Optimal Operation Model of Hydropower Reservoirs

The objective that is maximizing the total hydropower generation and complex constraint set from hydropower reservoirs constitute the optimal operation model of hydropower reservoirs.

2.1.1. Objectives

Please see equation bellow

$$\left\{\begin{array}{l}E=\max {\displaystyle \sum _{i=1}^N{\displaystyle \sum _{t=1}^T{P}_{i,t}}}\mathsf{\Delta }t\\ P_{i,t}=k_i{H}_{i,t}{Q}_{i,t}\end{array}\right.$$

(1)

where E is the total hydropower generation produced during operation periods; N is the number of hydropower reservoirs; T is the total operation periods; P_i,t, k_i,t, H_i,t, and Q_i,t are power output, efficiency coefficient, water head, and turbine discharge of reservoir i at period t, respectively; $\mathsf{\Delta }t$ is the operation interval.

2.1.2. Constraints

(1) Water balance Equations

$$\left\{\begin{array}{l}{V}_{i,t+1}=V_{i,t}+(I_{i,t}-Q_{i,t}^{total})\\ I_{i,t}=q_{i,t}+{\displaystyle \sum _{j=1}^{N_i}{O}_{j,t}}\\ O_{i,t}=Q_{i,t}^s+Q_{i,t}\end{array}\right.$$

(2)

where V_i,t, I_i,t, $Q_{i,t}^{total}$, q_i,t, and $Q_{i,t}^s$ are the storage volume, total inflow, total outflow, local inflow, and water spillage of reservoir i at period t, respectively; N_i is the number of upstream reservoirs directly connected to the reservoir i; O_j,t is the outflow from upstream reservoir j directly connected to the reservoir i at period t.

(2) Water head Equations

$$H_{i,t}=\frac{1}{2}(Z_{i,t}+Z_{i,t-1})-Z_{i,t}^d$$

(3)

where Z_i,t and $Z_{i,t}^d$ are the forebay water level and downstream water level of reservoir i at period t, respectively.

(3) Forebay water level limits

$$Z_{i,t}^{\min }\le Z_{i,t}\le Z_{i,t}^{\max }$$

(4)

where $Z_{i,t}^{\min }$ and $Z_{i,t}^{\max }$ are the minimum and maximum forebay water levels of reservoir i at period t, respectively.

(4) Forebay water level fluctuation limits

$$\left|Z_{i,t}-Z_{i,t-1}\right|\le \mathsf{\Delta }{Z}_i^{\max }$$

(5)

where $\mathsf{\Delta }{Z}_i^{\max }$ is the acceptable maximum fluctuation of forebay water level of reservoir i allowed for a timestep.

(5) Turbine discharge limits

$$Q_{i,t}^{\min }\le Q_{i,t}\le \min (Q_{i,t}^{\max },Q_{i,t}^c)$$

(6)

where $Q_{i,t}^{\min }$, $Q_{i,t}^{\max }$, and $Q_{i,t}^c$ are the minimum and maximum turbine discharges as well as reservoir discharge capacity of reservoir i at period t, respectively.

(6) Total discharge limits

$$Q_{i,t}^{total,\min }\le Q_{i,t}^{total}\le Q_{i,t}^{total,\max }$$

(7)

where $Q_{i,t}^{total,\min }$ and $Q_{i,t}^{total,\max }$ are the minimum and maximum total discharges of reservoir i at period t, respectively.

(7) Power output limits

$$P_{i,t}^{\min }\le P_{i,t}\le \min (P_{i,t}^{\max },P_{i,t}^e)$$

(8)

where $P_{i,t}^{\min }$, $P_{i,t}^{\max }$, and $P_{i,t}^e$ are the minimum and maximum power outputs as well as expected power outputs of reservoir i at period t, respectively.

(8) Total power output limits

$$\left\{\begin{array}{l}E=\max {\displaystyle \sum _{i=1}^N{\displaystyle \sum _{t=1}^T{P}_{i,t}}}\mathsf{\Delta }t\\ P_{i,t}=k_i{H}_{i,t}{Q}_{i,t}\end{array}\right.$$

(9)

where $P_t^{total,\hspace{0.33em}\min }$ and $P_t^{total,\hspace{0.33em}\max }$ are the minimum and maximum total power outputs of CHRs at period t.

(9) Initial and target forebay water level limits

$$\left\{\begin{array}{l}{Z}_{i,0}=Z_i^{initial}\\ Z_{i,T}=Z_i^{target}\end{array}\right.$$

(10)

where $Z_i^{initial}$ and $Z_i^{target}$ are the preset initial and target forebay water levels of reservoir i at period zero and T, respectively.

(10) Non-linear characteristic curves limits

$$\left\{\begin{array}{l}{V}_{i,t}=f_{i,1}(Z_{i,t})\\ Q_{i,t}^c=f_{i,2}(Z_{i,t-1})\\ Z_{i,t}^d=f_{i,3}(Q_{i,t}^{total})\\ P_{i,t}^e=f_{i,4}(H_{i,t},Q_{i,t})\\ Z_i^{target}=f_{i,1}^{-1}(V_i^s+V_{i,0})\end{array}\right.$$

(11)

where $f_{i,1}(\cdot )$, $f_{i,2}(\cdot )$, $f_{i,3}(\cdot )$, and $f_{i,4}(\cdot )$ are the nonlinear stage–storage curve, stage–discharge capacity curve, stage–downstream water level curve, and stage–discharge-head–output curve of reservoir i, respectively; $V_i^s$ is the required storage capacity of reservoir i during operation periods; $V_{i,0}$ is the storage volume of reservoir i at period zero.

(11) Storage and release water limits

$$V^s={\displaystyle \sum _{i=1}^N{V}_i^s}$$

(12)

where V^s is the total required storage capacity of multiple reservoir system.

It is noteworthy that the hybrid optimization algorithm combined dynamic programming (DP), successive approximation DP (DPSA) and progressive optimality algorithm (POA), called DP-POA-DPSA deterministic optimization method is used to solve the established optimal operation model of hydropower reservoirs, where penalty function method is used to deal with constraint violations.

2.2. Summary of Henry Gas Solubility Optimization (HGSO)

Henry’s law is a famous gas law in the field of physical chemistry formulated by famous chemical scientist Henry in 19th century [37]. Briefly, the law states that the amount of dissolved gas in a liquid is proportional to its partial pressure above the liquid when the gas and liquid reach equilibrium in the case of certain temperature. The key concept to understand Henry’s law is solubility, which changes with variation in temperature and pressure. Specifically, the increased temperature leads to the lower solubility of gas, conversely, increased pressure results in higher solubility of gas. By controlling the two factors of temperature and pressure, the highest equilibrium state with best gas can be achieved based on Henry’s law. Inspired by this phenomenon, Henry Gas Solubility Optimization (HGSO) algorithm was developed to mimic the physical and chemical process described by Henry’s law [38]. HGSO is a novel physics-inspired metaheuristic optimization algorithm based on population proposed in recent years and empirical study reveals that HGSO has a significant merit to balance between exploration and exploitation.

Similar to the majority of the well-known metaheuristic optimization algorithms, the optimization process of HGSO also includes two phases: the exploration phase and the exploitation phase. To achieve the balance between exploration and exploitation, the optimization process can be further divided into four stages: initialization stage, evaluation stage, updating stage, and termination stage. HGSO is no exception. In the initialization stage, the gas population are randomly generated with the population size N, lower bound vector lb, and upper bound vector ub as input, which can be formulated as follows:

$$g_i=lb+r_{1,i}(ub-lb)i=1, 2, \dots , N$$

(13)

where $g_i$ denotes the gas i; $r_{1,i}$ denotes the random number distributed in range from 0 to 1. It is noteworthy that the gas population are divided into several groups, denoted by l, with same Henry’s constant value. Then, each gas of the groups in population is evaluated based on the objective function in evaluation stage. Further, several key parameters such as Henry’s constant H, partial pressure P, and the ratio of enthalpy of dissolution to gas constant are randomly generated to update Henry’s coefficient and solubility for next generation, which can be summarized as follows:

$$H_l(it+1)=H_l(it)\exp (-C_l(\frac{1}{T(it)}-\frac{1}{T^{\theta }}))$$

(14)

$$S_{i,l}(it)=K{H}_l(it+1)P_{i,l}(it)$$

(15)

where $H_l(it)=r_2{c}_1$, $P_{i,l}=r_3{c}_2$, and $C_l=r_4{c}_3$, in which $r_2$, $r_3$, and $r_4$ are all random numbers in range from 0 to 1, $c_1$, $c_2$, and $c_3$ are constant values, which are set to 0.05, 100, and 0.01, respectively; it denotes the number of current iteration; $T(it)$ is the temperature in current iteration time and $T(it)=\exp (\frac{-it}{MaxIt})$, in which $MaxIt$ denotes the maximum number of iterations. Based on these parameters in current iteration, the gas positions of each group are updated in updating stage according to the following formulation:

$$g_{i,l}(it+1)=g_{i,l}(it)+r_5\gamma F(G_l(it)-g_{i,j}(it))+r_6\alpha F(S_{i,l}(it)G_l(it)-g_{i,j}(it))$$

(16)

$$\gamma =\beta \exp (-\frac{F(it)+\epsilon }{F_{i,l}(it)+\epsilon })$$

(17)

where $g_{i,l}$ denotes the gas i in group l; $r_5$ and $r_6$ are all random numbers in range from 0 to 1; $\gamma$ denotes the interaction ability of gases in the same group; $\alpha$, $\beta$, and $\epsilon$ are optimization parameters and are set to 1, 1, and 0.5, respectively; F denotes the direction parameter used to control the search direction of gas; $G_{i,l}(it)$ and $F_{i,l}(it)$ denote the best gas and its fitness in group l; $F(it)$ denotes the best fitness value obtained so far. In order to enhance the exploration capability of HGSO, the inferior gases are reupdated in the same way as initializing the population as follows:

$$N^{\prime }=N(c_4+r_7(c_5-c_4))$$

(18)

$$g_{i,l}^{\prime }=lb+r_8(ub-lb)$$

(19)

where $N^{\prime }$ denotes the number of inferior gases; $r_7$ and $r_8$ are random number in range from 0 to 1; $c_4=0.1$ and $c_5=0.1$; $g_{i,l}^{\prime }$ denotes the updated gas corresponding to gas i in group l. Finally, the HGSO algorithm is stopped when the termination criterion is met.

2.3. The Proposed MVQIHGSO

2.3.1. Multi-verse Optimizer (MVO)

The main inspiration for multi-verse optimizer comes from multi-verse theory, which has aroused great interest in the field of science, philosophy, and theology. In 2016, multi-verse optimizer was developed by Mirjalili and Lewis based on three simple yet important concepts including white holes, black holes, and wormholes to mimic the physic nexus of multiple universes in multi-verse theory [39]. Specially, in order to formulate the expansion of the universe, the higher expansion rate of the universe, the higher the probability of the existence of white holes, conversely, the lower expansion rate of the universe, the higher the probability of forming black holes.

The two key components guiding the optimization direction of the algorithm are white/black hole tunnel mechanism and wormhole mechanism. The former states that special transmission path can be established to realize the exchange of objects between different universes, such that the objects form black hole in the universe with low expansion rate will be replaced by the objects in the universe connected from the other end of the tunnel. In this manner, the inferior universe has the chance to explore the parallel space formed by multiple universe to evolve into the one with higher expansion rate. The mathematical model of white/black hole tunnel can be formulated as follows:

$$u_{i,j}=\left\{\begin{array}{l}{u}_{k,j} r_1<normr(u_i)\\ u_{i,j} r_1\ge normr(u_i)\end{array}\right.$$

(20)

where $u_i$ denotes the universe i; $u_{i,j}$ denotes the object j of universe i, wherein $j=1,2,\dots ,D$ and D denotes the number of objects in a universe. $normr(\cdot )$ aims to normalize the expansion rate of universe to a length of 1; $r_1$ denotes the random number uniformly distributed on [0, 1]. It is noteworthy that universe and object correspond to the individual and decision variable in optimization algorithm, respectively.

Additionally, the wormhole belongs to the universe with highest expansion rate. In order to share the information of objects in wormhole with other universes, objects in wormholes have a high probability of being randomly transmitted to other universes. In this way, the universe diversity is enhanced, avoiding slow expansion or stopping expansion during the process of evolution. The wormhole mechanism can be mathematically described as follows:

$$u_{i,j}=\left\{\begin{array}{l}{U}_j+TDR(l{b}_j+r_2(u{b}_j-l{b}_j)) r_3<0.5 \mathrm{and} r_4<WEP\\ U_j-TDR(l{b}_j+r_2(u{b}_j-l{b}_j)) r_3\ge 0.5 \mathrm{and} r_4<WEP\\ u_{i,j} r_4\ge 0.5\end{array}\right.$$

(21)

where $U_j$ denotes the object j of the best universe; $r_2$, $r_3$, and $r_4$ denote the random number uniformly distributed on [0, 1], respectively; TDR and WEP denotes travelling distance rate and wormhole existence probability, respectively, wherein the former is used to control the transmission range of objects in the best universe and the latter is applied to reflect the possibility of the existence of wormhole. The detailed information of MVO can be referred to the original paper [39].

2.3.2. Quadratic Interpolation Strategy (QI)

Quadratic interpolation has been proved to be a commonly and effective method to improve the exploitation ability of swarm intelligence algorithm [40,41,42]. QI is essentially a one-dimensional optimization method, which can easily and quickly discover the optimal solution in the context of low-dimensional space. Its core idea is to construct a quadratic curve to simulate the objective function around several candidate solutions, so that the optimal function value can be calculated in a derivative-free way, which greatly reduces the cost of solving the local objective value. In other words, QI method tries to find better promising search agent in the vicinity of individuals evolved so far, leading to improve the population diversity and enhance the exploitation ability of the algorithm, so that it can increase the probability of guiding the search agent toward the global optimal region.

Generally, QI method is the last step before evolution into the next generation. In this step, three individuals will be selected to construct a quadratic curve to simulate the objective function. In this manner, the extreme point of the approximate objective function can be easily obtained by quadratic function. For example, $s_a=(s_{a,1},s_{a,2},\dots ,s_{a,j},\dots ,s_{a,D})$, $s_b=(s_{b,1},s_{b,2},\dots ,s_{b,j},\dots ,s_{b,D})$, and $s_c=(s_{c,1},s_{c,2},\dots ,s_{c,j},\dots ,s_{c,D})$ are the three individuals selected to form the quadratic curve, the extreme point $s_p=(s_{p,1},s_{p,2},\dots s_{p,j},\dots ,s_{p,D})$ of the quadratic, which is also the promising possible solution based on the three selected individuals, can be formulated as follows:

$$s_{p,j}=\frac{({(s_{a,j})}^2-{(s_{c,j})}^2)f(s_b)+({(s_{c,j})}^2-{(s_{b,j})}^2)f(s_a)+({(s_{b,j})}^2-{(s_{a,j})}^2)f(s_c)}{2((s_{a,j}-s_{c,j})f(s_b)+(s_{c,j}-s_{b,j})f(s_a)+(s_{b,j}-s_{a,j})f(s_c))}$$

(22)

where $j=1,2,\dots ,D$; $f(\cdot )$ denotes the fitness value of a specific individual. In practice, the individuals in current iteration are ranked from small to large based on their fitness values. Then, all three consecutive individuals are selected to perform the QI method. As a result, for the population with size N where $N\ge 3$, the QI method is implemented N-2 times. Based on the greedy law, if the resulting individual by QI method performs better than the selected three individuals, the first individual will be replaced.

As can be seen from the above analysis, QI can be regarded as a crossover operator to an extent. In other words, the potential information hidden in individuals is fully utilized to improve the performance of inferior solution, so that the promising solution with high quality can be discovered.

2.3.3. The Proposed Enhanced HGSO Algorithm

In the proposed MVQIHGSO algorithm, the MVO operator and QI strategy are MVO operator is applied to enhance the exploration capability and improve the balance between exploration and exploitation of basic HGSO. At the same time, QI strategy is used to find the promising candidate solution around the search space formed by the current individual. In other words, the inferior solution is replaced with the solution obtained by QI strategy to increase the diversity and quality of population, leading to faster convergence rate and higher exploitation accuracy. The pseudo-code of MVQIHGSO is given in Algorithm 1.

Algorithm 1. Detailed information of MVQIHGSO algorithm.

Pseudocode of the MVQIHGSO 01: Inputs: the population size N; the maximum number of iteration MaxIt; the constant values including c1, c2, c3 in Equations (2) and (3), c4 and c5 in Equation (6), α, β, and ε in Equatios (4) and (5); the number of groups l 02: Initialize the gas population within the lower and upper boundary 03: Divide the gas population into specific groups with the same Henry’s constant value 04: Evaluate the gas of each group in the population 05: Obtain the best gas of the whole population and the l best gases corresponding to l groups 06: for it from 1 to MaxIt do 07: Generate random number r within [0, 1] 08: if r smaller than 0.5 do 09:  Update the gas position by basic HGSO 10: else do 11:  Update the gas position by MVO 12: end if 13: Sort the updated gas population by HGSO and MVO 14: Updating the gas population from 1 to N-2 by QI strategy based on greedy law 15: Evaluate the fitness of the final updated gas population 16: Update the best gas position obtained so far and the l best gases corresponding to l groups 17: end for 18: Outputs: the best global gas position and the corresponding optimal fitness value

2.4. Support Vector Machine

Support vector machine (SVM) is a famous machine learning algorithm based on statistical learning theory and the principle of structural risk minimization that has received the most attention. SVM was firstly proposed by Cortes and Vapnik with excellent performance in term of high generalization ability [32]. Recently, there has been a growing interest from many scholars and researchers in applying SVM in various fields for classification and regression analysis such as streamflow forecasting, classification of image, text and hypertext categorization, classification of satellite data, and classification of proteins, which have achieved lots of advances. The main idea of SVM is transforming the nonlinear input space to a high-dimensional feature space to construct a hyper-plane by nonlinear mapping [43].

The operation rule derivation belongs to regression problem, so epsilon-insensitive SVM (ɛ-SVM) regression is applied in the present study. ɛ-SVM was designed by Drucker in 1997 to solve the regression problem [44]. Suppose the truth yet unknown function of input space x as G(x) and F(x, w) is a family of functions controlled by weight vector w, the aim of SVM is find a best w denoted as $\hat{w}$ to minimize a measure of error between G(x) and a specific approximate function of F(x, w). The optimal value of $\hat{w}$ depends on the primal loss function given as follows:

$$L=\left\{\begin{array}{l} 0 \mathrm{if} \left|y_i-F(x_i,\widehat{w})\right|<\epsilon \\ \left|y_i-F(x_i,\widehat{w})\right|-\epsilon \mathrm{otherwise}\end{array}\right.$$

(23)

where x_i is the input vector of sample point i; y_i is the observation value corresponding to sample point i; $F(x_i,\widehat{w})$ denotes the predicted value of sample point i; $\epsilon$ is defined as a cube so that the loss is zero in the case of $F(x_i,\widehat{w})$ is within the tube, while the loss is the magnitude of the difference between $F(x_i,\widehat{w})$ and the radius $\epsilon$ of the tube [44]. In order to solve the optimization problem in term of the above primal loss function, the slack variables ${\xi }^{\ast }$ and $\xi$ are introduced to minimize:

$$U⟮{\displaystyle \sum _{i=1}^N{\xi }_i^{\ast }+{\displaystyle \sum _{i=1}^N{\xi }_i}}⟯+\frac{1}{2}(w^{t}w)$$

(24)

where $U(\cdot )$ is the function of ${\xi }^{\ast }$ and $\xi$, which emphasizes more on the error when U is large, and vice versa, it emphasizes more on the norm of the weights resulting in a better generalization; N is the number of training vectors. The corresponding constraints are as follows:

$$\begin{array}{l}{y}_i-(w^t{v}_i)-b=\epsilon +{\xi }_i\\ (w^t{v}_i)+b-y_i\le \epsilon +{\xi }_i^{\ast }\\ {\xi }_i^{\ast }\ge 0\\ {\xi }_i\ge 0\end{array}$$

(25)

where $v_i$ is the training vector i; b denotes a group of bias. On this basis, the Lagrangian can be formulated based on the Lagrange multipliers ${\gamma }_i$ and ${\alpha }_i$ as follows:

$$\begin{array}{l}L=\frac{1}{2}(w^{t}w)+U⟮{\displaystyle \sum _{i=1}^N{\xi }_i^{\ast }+{\displaystyle \sum _{i=1}^N{\xi }_i}}⟯-{\displaystyle \sum _{i=1}^N{\alpha }_i^{\ast }[y_i-(w^t{v}_i)-b+\epsilon +{\xi }_i^{\ast }]}\\ -{\displaystyle \sum _{i=1}^N{\alpha }_i^{\ast }[(w^t{v}_i)+b-y_i+\epsilon +{\xi }_i]}-{\displaystyle \sum _{i=1}^N({\gamma }_i^{\ast }{\xi }_i^{\ast }+{\gamma }_i{\xi }_i)}\end{array}$$

(26)

As can be seen from Equation (14) that in addition to parameters w, b, and ${\xi }_i$, Lagrangian is also related to the Lagrange multipliers, which makes the direct solution of the parameters of ɛ-SVM with respect to w, b, and ${\xi }_i$ by taking the partial derivative of the Lagrangian impossible. Fortunately, the saddle point of L by differentiating with respect to w, b, and ${\xi }_i$ leads to the equivalent maximization of the dual space, which can be formulated as follows:

$$W(\alpha ,{\alpha }^{\ast })=-\epsilon {\displaystyle \sum _{i=1}^N({\alpha }_i^{\ast }+{\alpha }_i})+{\displaystyle \sum _{i=1}^N{y}_i({\alpha }_i^{\ast }-{\alpha }_i})-\frac{1}{2}{\displaystyle \sum _{i,j=1}^N({\alpha }_i^{\ast }-{\alpha }_i})({\alpha }_j^{\ast }-{\alpha }_j){(v^{t}v+1)}^p$$

(27)

with the constraints:

$$\begin{array}{l}0\le {\alpha }_i\le U\hspace{1em}0\le {\alpha }_i^{\ast }\le Ui=1, \dots , N\\ {\displaystyle \sum _{i=1}^N{\alpha }_i^{\ast }}={\displaystyle \sum _{i=1}^N{\alpha }_i}\end{array}$$

(28)

As a result, the prediction of $y^{(p)}$ corresponding to new input vector $x^{(p)}$ can be formulated as follows:

$$y^{(p)}={\displaystyle \sum _{i=1}^N({\alpha }_i^{\ast }-{\alpha }_i)(v_i^t}{x}^{(p)}+1{)}^p+b$$

(29)

In order to effectively and efficiently solve various linear and nonlinear regression problems, the most critical component of SVM, namely kernel function, is described here. The application of kernel function ensures that the vector dot product of high-dimensional space can be directly calculated in low-dimensional space according to the Hilbert–Schmidt theory [45], greatly reducing calculation complexity and avoiding the curse of dimensional. In this manner, the performance of SVM either in classification and regression is significantly improved. The operation of kernel function is as follows:

$$K(x_1,x_2)=〈\varphi (x_1),\varphi (x_2)〉$$

(30)

where $K(\cdot )$ is the kernel function; $\varphi (x)$ is a transformation that maps $x$ to a high-dimensional space; $〈\cdot 〉$ denotes dot product operator. The commonly used kernel functions in SVM include linear function, polynomial function, radial basis function (RBF), and sigmoidal function. Generally, RBF can achieve significant superiority against other kernel functions, which can be expressed as follows:

$$K_{RBF}(x_1,x_2)=\exp (-\gamma {‖x_1-x_2‖}^2) \gamma > 0$$

(31)

where $\gamma$ is one hyperparameter of SVM.

It is noteworthy that data processing is needed to standardize the input samples to eliminate the numerical calculation instability caused by feature dimension and improve the convergence speed. The standardization of input sample x can be expressed as follows:

$$x^{\ast }=\frac{x-u}{\sigma }$$

(32)

where $x^{\ast }$ is the standardized sample; $u$ and $\sigma$ are the mean vector and standard deviation vector, respectively.

In summarize, the MVQIHGSO and SVM are combined to construct a novel optimal operation rule derivation method, which is called MVQIHGSO-SVM model. The proposed MVQIHGSO is used to optimize the hyperparameters of SVM, the flowchart of MVQIHGSO-SVM is shown in Figure 1. Furthermore, several commonly used validation indices are used to comprehensively evaluate the performance of the proposed MVQIHGSO-SVM model. The validation indices are composed of coefficient of determination (R²), root mean square error (RMSE), mean absolute error (MAE) and mean absolute percentage error (MAPE). Different indices can validate the performances of the model in different dimensions. Specifically, R² is used to characterize the linear correlation between the observed data and its predicted data; RMSE is adopted to indicate the average degree of dispersion between the observed data and its predicted data; MAE and MAPE are applied to quantify the deviation of the predicted data from the observed data. The formulation of validation indices are as follows:

$$R^2=1-\frac{{\displaystyle \sum _{i=1}^n{(y_i-{\widehat{y}}_i)}^2}}{{\displaystyle \sum _{i=1}^n{(y_i-{\overline{y}}_i)}^2}}$$

(33)

$$RMSE=\sqrt{\frac{1}{n}{\displaystyle \sum _{i=1}^n{(y_i-{\widehat{y}}_i)}^2}}$$

(34)

$$MAE=\frac{1}{n}{\displaystyle \sum _{i=1}^n\left|y_i-{\widehat{y}}_i\right|}$$

(35)

$$MAPE=\frac{100}{n}{\displaystyle \sum _{i=1}^n\left|\frac{y_i-{\widehat{y}}_i}{y_i}\right|}$$

(36)

It is noteworthy that the smaller these indicators, the better the performance of the model.

3. Experimental Evaluation and Discussion

In order to verify the performance of the proposed MVQIHGSO algorithm, a comprehensive set of benchmark functions including 23 widely used functions in the field of intelligence algorithm are adopted as test bed. These benchmark function set includes three types: unimodal (F1–F7), multimodal (F8–F13), and fixed-dimension multimodal (F14–F23), which are all minimization optimization problems. F1–F7 are used to test the convergence speed and exploitation accuracy of algorithms because of only one local or global optima exists in the search space; F8–F13 are used to test the exploration capacity when facing many local optimal solutions; F14–F23 are used to test the ability of local optima avoidance and the balance between exploration and exploitation of the algorithms. More detailed information about these benchmark functions including function expressions, dim and boundary of decision variables, and the theoretical global optima are presented in

Table 1. Description of unimodal benchmark functions.

Function	Dim	Range	fmin
$f_1(x)={\displaystyle {\sum }_{i=1}^n{x}_i^2}$	30	[−100, 100]	0
$f_2(x)={\displaystyle {\sum }_{i=1}^n\left|x_i^{}\right|}+{\prod }_{i=1}^n\left|x_i^{}\right|$	30	[−10, 10]	0
$f_3(x)={{\displaystyle {\sum }_{i=1}^n\left({\displaystyle {\sum }_{j-1}^i{x}_j}\right)}}^2$	30	[−100, 100]	0
$f_4(x)={\max }_i\left\{\left|x_i\right|,1\le i\le n\right\}$	30	[−100, 100]	0
$f_5(x)={\displaystyle {\sum }_{i=1}^{n-1}\left[100{\left(x_{i+1}-x_i^2\right)}^2+{\left(x_i-1\right)}^2\right]}$	30	[−30, 30]	0
$f_6(x)={\displaystyle {\sum }_{i=1}^n{\left(\left[x_i+0.5\right]\right)}^2}$	30	[−100, 100]	0
$f_7(x)={\displaystyle {\sum }_{i=1}^{n}i{x}_i^4}+random\left[0,1\right)$	30	[−1.28, 1.28]	0

,

Table 2. Description of multimodal benchmark functions.

Function	Dim	Range	fmin
$f_8(x)={\displaystyle {\sum }_{i=1}^n-x_i^{}}\sin (\sqrt{\left|x_i\right|})$	30	[−500, 500]	−418.9829 × 5
$f_9(x)={\displaystyle {\sum }_{i=1}^n\left[x_i^2-10\cos (2\pi x_i)+10\right]}$	30	[−5.12, 5.12]	0
$f_{10}(x)=-20\exp ⟮-0.2\sqrt{\frac{1}{n}{\displaystyle {\sum }_{i=1}^n{x}_i^2}}⟯-\exp ⟮\frac{1}{n}{\displaystyle {\sum }_{i=1}^n\cos (2\pi x_i)}⟯+20+e$	30	[−32, 32]	0
$f_{11}(x)=\frac{1}{4000}{\displaystyle {\sum }_{i=1}^n{x}_i^2}-{\prod }_{i=1}^n\cos ⟮\frac{x_i}{\sqrt{i}}⟯+1$	30	[−600, 600]	0
$\begin{array}{l}{f}_{12}(x)=\frac{\pi }{n}\left\{10\sin (\pi y_1)+{\displaystyle {\sum }_{i=1}^{n-1}{(y_i-1)}^2\left[1+10{\sin }^2(\pi y_{i+1})\right]+{(y_n-1)}^2}\right\}+{\displaystyle {\sum }_{i=1}^{n}u(x_i,10,100,4)}\\ y_i=1+\frac{x_i+1}{4}\\ u(x_i,a,k,m)=\left\{\begin{array}{l}k{\left(x_i-a\right)}^m x_i>a\\ 0 -a<x_i<a\\ k{\left(-x_i-a\right)}^m x_i<-a\end{array}\right.\end{array}$	30	[−50, 50]	0
$f_{13}(x)=0.1\left\{{\sin }^2(3\pi x)+{\displaystyle {\sum }_{i=1}^n{(x_i-1)}^2\left[1+{\sin }^2(3\pi x_i+1)\right]+{(x_n-1)}^2\left[1+{\sin }^2(2\pi x_n)\right]}\right\}+{\displaystyle {\sum }_{i=1}^{n}u(x_i,5,100,4)}$	30	[−50, 50]	0

and

Table 3. Description of fixed-dimension multimodal benchmark functions.

Function	Dim	Range	fmin
$f_{14}(x)={⟮\frac{1}{500}+{\displaystyle {\sum }_{j=1}^{25}\frac{1}{j+{\displaystyle {\sum }_{i=1}^2{(x_i-a_{ij})}^6}}}⟯}^{-1}$	2	[−65, 65]	1
$f_{15}(x)={{\displaystyle {\sum }_{i=1}^{11}\left[a_i-\frac{x_1(b_i^2+b_i{x}_2)}{b_i^2+b_i{x}_3+x_4}\right]}}^2$	4	[−5, 5]	0.00030
$f_{16}(x)=4x_1^2-2.1x_1^4+\frac{1}{3}{x}_1^6+x_1{x}_2-4x_2^2+4x_2^4$	2	[−5, 5]	−1.0316
$f_{17}(x)=⟮x_2-\frac{5.1}{4{\pi }^2}{x}_1^2+\frac{5}{\pi }{x}_1-6⟯+10⟮1-\frac{1}{8\pi }⟯\cos x_1+10$	2	[−5, 5]	0.398
$f_{18}(x)=\left[1+{(x_1+x_2+1)}^2(19-14x_1+3x_1^2-14x_2+6x_1{x}_2+3x_2^2)\right]\times \left[30+{(2x_1-3x_2)}^2\times (18-32x_1+12x_1^2+48x_2-36x_1{x}_2+27x_2^2)\right]$	2	[−2, 2]	3
$f_{19}(x)=-{\displaystyle {\sum }_{i=1}^4{c}_i}\exp \left(-{\displaystyle {\sum }_{j=1}^3{a}_{ij}{(x_j-p_{ij})}^2}\right)$	3	[1, 3]	−3.86
$f_{20}(x)=-{\displaystyle {\sum }_{i=1}^4{c}_i}\exp \left(-{\displaystyle {\sum }_{j=1}^6{a}_{ij}{(x_j-p_{ij})}^2}\right)$	6	[0, 1]	−3.32
$f_{21}(x)=-{{\displaystyle {\sum }_{i=1}^5\left[(X-a_i){(X-a_i)}^T+c_i\right]}}^{-1}$	4	[0, 10]	−10.1532
$f_{22}(x)=-{{\displaystyle {\sum }_{i=1}^7\left[(X-a_i){(X-a_i)}^T+c_i\right]}}^{-1}$	4	[0, 10]	−10.4028
$f_{23}(x)=-{{\displaystyle {\sum }_{i=1}^{10}\left[(X-a_i){(X-a_i)}^T+c_i\right]}}^{-1}$	4	[0, 10]	−10.5363

.

In addition, for a general and fair comparison, the proposed MVQIHGSO are compared with the basic HGSO and several state-of-art metaheuristic optimization algorithms including particle swarm optimization (PSO) [46], differential evolution algorithm (DE) [47], multi-verse optimizer (MVO) [36], sine cosine algorithm (SCA) [48], opposition-based sine cosine algorithm (OBSCA) [49], grey wolf optimizer (GWO) [50], and improved grey wolf optimizer (IGWO) [51]. It is noteworthy that the parameters settings of the comparison algorithms are the same as those used in the original literature. At the same time, without the loss of generality, for the specific common parameters, for example, the population size, the dim of decision variables, and the maximum number of iterations are set to 50, 30, and 1000. Specially, in order to reduce the sudden impact of the randomness of the algorithm on the overall performance, all the optimization algorithms are carried out with 30 runs independently for each benchmark function to analyze the statistical performance of the algorithms.

3.1. Statistical Results and Analysis

In this subsection, the aforementioned algorithms were run 30 times on each benchmark function, such that the statistical results including average and standard deviation were calculated based on the optimization results. The corresponding statistical results were recorded and presented in

Table 4. Statistical results of the proposed MVQIHGSO and other state-of-art algorithms.

Functions		Indicator	MVQIHGSO	HGSO	PSO	DE	MVO	SCA	OBSCA	GWO	IGWO
Unimodal	F1	Mean	0.00	1.71 × 10−71	2.33 × 10−40	5.29 × 10−12	1.79 × 10−1	2.74 × 10−3	1.11 × 10−27	2.48 × 10−70	2.57 × 10−71
		Std	9.34 × 10−71	0.00	1.04 × 10−39	2.22 × 10−12	5.80 × 10−2	7.88 × 10−3	4.22 × 10−27	7.10 × 10−70	6.07 × 10−71
	F2	Mean	4.10 × 10−185	1.25 × 10−43	8.79 × 10−3	4.65 × 10−8	3.55	3.17 × 10−6	4.21 × 10−25	6.67 × 10−41	1.95 × 10−42
		Std	6.83 × 10−43	0.00	3.15 × 10−2	1.23 × 10−8	1.81 × 101	4.67 × 10−6	2.18 × 10−24	7.20 × 10−41	3.62 × 10−42
	F3	Mean	0.00	1.04 × 10−69	1.09 × 10−1	2.41 × 104	1.83 × 101	2.10 × 103	1.32 × 10−3	7.80 × 10−19	8.61 × 10−13
		Std	5.66 × 10−69	0.00	1.28 × 10−1	2.77 × 103	6.90	1.73 × 103	5.17 × 10−3	3.92 × 10−18	4.38 × 10−12
	F4	Mean	2.06 × 10−72	1.51 × 10−183	1.09 × 10−1	2.00	6.33 × 10−1	1.36 × 101	1.55 × 10−5	1.15 × 10−17	7.93 × 10−15
		Std	1.13 × 10−71	0.00	7.10 × 10−2	2.44 × 10−1	2.77 × 10−1	9.40	2.88 × 10−5	1.69 × 10−17	6.16 × 10−15
	F5	Mean	2.85 × 101	2.80 × 101	4.06 × 101	4.65 × 101	3.08 × 102	5.54 × 102	2.79 × 101	2.65 × 101	2.25 × 101
		Std	2.70 × 10−1	6.13 × 10−1	2.77 × 101	2.31 × 101	5.98 × 102	2.18 × 103	3.06 × 10−1	7.80 × 10−1	2.79 × 10−1
	F6	Mean	1.74 × 10−1	3.44	3.63 × 10−23	4.92 × 10−12	1.66 × 10−1	4.29	4.10	4.13 × 10−1	1.01 × 10−5
		Std	6.39 × 10−2	5.26 × 10−1	1.84 × 10−22	1.84 × 10−12	4.89 × 10−2	3.92 × 10−1	2.73 × 10−1	2.57 × 10−1	2.39 × 10−6
	F7	Mean	6.76 × 10−5	7.92 × 10−4	9.32 × 10−3	2.59 × 10−2	1.22 × 10−2	2.35 × 10−2	1.47 × 10−3	4.72 × 10−4	8.64 × 10−4
		Std	4.35 × 10−4	4.84 × 10−5	4.00 × 10−3	4.70 × 10−3	5.04 × 10−3	2.54 × 10−2	1.03 × 10−3	3.15 × 10−4	3.80 × 10−4
Multimodal	F8	Mean	−1.02 × 104	−2.64 × 105	−6.68 × 103	−1.25 × 104	−8.18 × 103	−3.97 × 103	−4.08 × 103	−6.09 × 103	−9.62 × 103
		Std	1.08 × 103	6.38 × 105	5.86 × 102	8.39 × 101	7.31 × 102	2.69 × 102	2.26 × 102	7.63 × 102	1.29 × 103
	F9	Mean	0.00	0.00	4.67 × 101	6.20 × 101	1.05 × 102	1.29 × 101	0.00	1.78 × 10−1	1.39 × 101
		Std	0.00	0.00	1.45 × 101	5.96	3.31 × 101	1.98 × 101	0.00	8.08 × 10−1	6.96
	F10	Mean	1.01 × 10−15	1.72 × 10−15	6.36 × 10−1	6.01 × 10−7	7.92 × 10−1	1.11 × 101	1.09 × 10−1	1.26 × 10−14	9.06 × 10−15
		Std	1.53 × 10−15	6.49 × 10−16	7.67 × 10−1	1.10 × 10−7	7.47 × 10−1	9.72	5.12 × 10−1	2.97 × 10−15	2.31 × 10−15
	F11	Mean	0.00	0.00	1.66 × 10−2	7.82 × 10−11	4.49 × 10−1	1.72 × 10−1	4.67 × 10−11	4.53 × 10−4	1.89 × 10−3
		Std	0.00	0.00	2.13 × 10−2	1.51 × 10−10	8.69 × 10−2	2.35 × 10−1	2.56 × 10−10	2.48 × 10−3	4.56 × 10−3
	F12	Mean	7.73 × 10−4	3.43 × 10−1	4.49 × 10−2	6.66 × 10−13	8.71 × 10−1	1.17	4.48 × 10−1	3.01 × 10−2	7.46 × 10−7
		Std	3.65 × 10−4	1.18 × 10−1	7.04 × 10−2	3.91 × 10−13	8.25 × 10−1	1.89	9.21 × 10−2	2.30 × 10−2	2.44 × 10−7
	F13	Mean	2.00 × 10−2	2.49	2.07 × 10−2	3.01 × 10−12	3.78 × 10−2	3.20	2.28	3.05 × 10−1	1.63 × 10−2
		Std	8.67 × 10−3	3.23 × 10−1	3.65 × 10−2	1.72 × 10−12	1.87 × 10−2	1.52	1.47 × 10−1	2.03 × 10−1	3.70 × 10−2
Fixed-dimension multimodal	F14	Mean	9.98 × 10−1	1.14	2.58	9.98 × 10−1	9.98 × 10−1	1.33	1.20	3.22	9.98 × 10−1
		Std	1.61 × 10−12	2.87 × 10−1	2.01	0.00	6.12 × 10−12	7.52 × 10−1	6.05 × 10−1	3.54	4.12 × 10−17
	F15	Mean	3.08 × 10−4	3.53 × 10−4	3.84 × 10−4	6.49 × 10−4	3.33 × 10−3	9.72 × 10−4	7.45 × 10−4	4.35 × 10−3	3.34 × 10−4
		Std	4.16 × 10−8	4.13 × 10−5	2.95 × 10−4	9.37 × 10−5	6.80 × 10−3	4.47 × 10−4	1.14 × 10−4	8.14 × 10−3	1.43 × 10−4
	F16	Mean	−1.03	−1.03	−1.03	−1.03	−1.03	−1.03	-1.03	−1.03	−1.03
		Std	5.34 × 10−12	2.10 × 10−5	6.78 × 10−16	6.78 × 10−16	4.39 × 10−8	1.50 × 10−5	7.18 × 10−7	2.78 × 10−9	6.78 × 10−16
	F17	Mean	0.398	0.399	0.398	0.398	0.398	0.399	0.398	0.398	0.398
		Std	1.49 × 10−10	9.22 × 10−4	0.00	0.00	7.35 × 10−8	6.37 × 10−4	2.89 × 10−4	3.25 × 10−7	0.00
	F18	Mean	3.00	3.00	3.00	3.00	3.00	3.00	3.00	3.00	3.00
		Std	2.86 × 10−10	1.13 × 10−5	1.50 × 10−15	1.90 × 10−15	6.23 × 10−7	1.51 × 10−5	4.30 × 10−6	3.90 × 10−6	7.14 × 10−16
	F19	Mean	−3.86	−3.86	−3.86	−3.86	−3.86	−3.86	−3.86	−3.86	−3.86
		Std	1.53 × 10−9	2.17 × 10−3	2.71 × 10−15	2.71 × 10−15	2.67 × 10−7	2.89 × 10−3	1.76 × 10−3	1.30 × 10−3	2.71 × 10−15
	F20	Mean	−3.32	−3.12	−3.29	−3.31	−3.29	−2.96	−3.16	−3.29	−3.31
		Std	3.02 × 10−2	7.03 × 10−2	5.54 × 10−2	1.38 × 10−15	5.56 × 10−2	3.10 × 10−1	4.12 × 10−2	5.20 × 10−2	3.02 × 10−2
	F21	Mean	−10.2	−4.89	−6.15	−9.96	−7.28	−4.28	−9.30	−9.31	−9.82
		Std	1.93	8.91 × 10−2	3.44	1.68 × 10−3	2.79	2.16	1.07 × 10−1	1.92	1.28
	F22	Mean	−10.4	−4.91	−9.06	−10.0	−8.58	−4.05	−10.2	−10.2	−10.4
		Std	1.35	1.15 × 10−1	2.77	5.07 × 10−8	2.90	2.26	1.44 × 10−1	9.63 × 10−1	7.32 × 10−9
	F23	Mean	−10.5	−4.96	−8.11	−10.5	−9.56	−4.97	−10.4	−10.4	−10.5
		Std	2.59 × 10−5	1.02 × 10−1	3.55	1.49 × 10−13	2.25	1.70	1.04 × 10−1	9.79 × 10−1	1.09 × 10−14

, where the best statistical results with respect to average values of algorithm are presented in bold. As shown in Table 4 that MVQIHGSO is able to provide very competitive results and has the smallest average values on 16 out of the 23 benchmark functions, more than half of them. More specifically, the analysis of specific performance is as follows.

For unimodal functions, the proposed algorithm has achieved significant superiority against other well-known algorithms in F1, F2, F3, and F7, indicating the superior performance of MVQIHGSO in term of finding the global optimum in local search space during exploitation phase. This is because the quadratic interpolation method is imbedded into the basic HGSO algorithm, such that the ability of the basic HGSO to exploit the optimum is enhanced. For multimodal functions, the proposed algorithm also outperforms other algorithms in F9, F10, and F11, indicating its good exploration capacity in the case of a large number of local optima because of the introduced white/black hole mechanism from MVO algorithm to solve local optima stagnation by sudden change. For fixed-dimension multimodal functions, the proposed algorithm achieves the best performance in all the benchmark functions compared with other algorithms caused by the incorporated wormhole mechanism with adaptive controlling parameters to balance exploration and exploitation. The above analysis based on statistical results reveal that the resulting algorithm has merits with respect to robust and statistically sound in term of exploration phase and exploitation phase.

3.2. Non-Parameter Test Results and Analysis

Friedman test is a non-parametric statistical test developed by Milton Friedman, which is often used to detect differences in treatments across multiple test attempts [52,53]. For the performance comparison of metaheuristic algorithms, it uses rank information based on the statistical results of several runs to examine significant differences in multiple population distributions [54]. Therefore, in this subsection, Friedman test was first used to benchmark the overall performance of MVQIHGSO and other state-of-art algorithms. In this manner, the average ranks of the algorithms are obtained in

Table 5. Statistical results of the proposed MVQIHGSO and other state-of-art algorithms from Friedman test.

Algorithms	Friedman Ranks	Final Ranks
MVQIHGSO	2.543	1
HGSO	4.608	3
PSO	5.652	7
DE	4.608	4
MVO	6.456	8
SCA	7.804	9
OBSCA	5.260	6
GWO	4.695	5
IGWO	3.369	2

. As shown in Table 5 that the Friedman rank of the proposed MVQIHGSO is 2.543, which is the smallest value among all the algorithms. The second-best algorithm is IGWO with the Friedman rank 3.369. These results show that the overall performance of the proposed algorithm is better in terms of exploration and exploitation.

Additionally, in order to further explore the performance difference in each run to exclude the case that superiority and dissimilarity occurs by chance, another non-parameter test named Wilcoxon singed-rank test [55] was also carried out to determine whether there is a statistical significant difference between the proposed MVQIHGSO and other state-of-art algorithms in this subsection. The null hypothesis is that the differences have a median of zero and the alternative hypothesis that the median is not equal to zero for a two-sided test or greater (or smaller) than zero for a one-sided test. The results of Wilcoxon singed-rank test with a significance level of 5% is represented in

Table 6. Statistical results of the proposed MVQIHGSO and other state-of-art algorithms from Wilcoxon test ($p\ge 0.05$).

Compared Algorithms	Unimodal Functions	Multimodal Functions	Fixed-Dimension Functions
MVQIHGSO vs. HGSO	2.5940 × 10−8	0.1012	7.4567 × 10−4
MVQIHGSO vs. PSO	1.4838 × 10−12	4.1333 × 10−10	0.8573
MVQIHGSO vs. DE	6.6342 × 10−19	2.5731 × 10−4	0.0916
MVQIHGSO vs. MVO	1.7618 × 10−32	5.2700 × 10−31	0.0273
MVQIHGSO vs. SCA	9.6394 × 10−27	6.4376 × 10−25	3.2725 × 10−5
MVQIHGSO vs. OBSCA	1.0240 × 10−9	4.1657 × 10−4	0.0227
MVQIHGSO vs. GWO	5.5859 × 10−7	0.0072	0.0654
MVQIHGSO vs. IGWO	7.5243 × 10−5	0.0011	0.0402

. As shown in Table 6, the proposed algorithm performs significantly better than other algorithms in the case where the p value is less than 0.05. This can strongly demonstrate the potential capabilities of the proposed algorithm in coping with the complex optimization problems in real world.

From the analysis of Friedman test and Wilcoxon rank sum test, it can be concluded that the MVQIHGSO is able to avoid local optima stagnation and find the global optimal solution more accurately. This is due to the introduction of MVO operator and QI strategy for improving the position updating mechanism of the basic HGSO. In this way, the convergence speed can also be improved, which can be found in Figure 2. Figure 2 shows that the convergence speed of MVQIHGSO is significantly higher than other algorithms.

4. Case Study

4.1. Study Region

The main stream of the Jinsha River is the largest hydropower plant planned and constructed in China. The lower reaches of the Jinsha River (from the Yalong River Estuary to Yibin) have the richest hydropower resources, with a length of 782 km and a drop of 729 m. Xiluodu and Xiangjiaba (XLD and XJB) hydropower reservoirs, located in the lower reaches of the Jinsha River, are China’s third and fifth largest hydropower stations with a total installed capacity of 18,600 MW, equivalent to building a new Three Gorges hydropower station. The landscape of XLD and XJB hydropower stations are shown in Figure 3, and the characteristic parameters are given in

Table 7. Characteristics information of XLD-XJB cascade reservoirs.

Characteristics	Hydropower Stations		Units
	Xiluodu	Xiangjiaba
Completion date	2013	2012	-
Watershed area	0.45	0.45	million km2
Dead water level	540	370	m
Flood control limited water level	560	370	m
Normal water level	600	380	m
Regulated storage	6.46	0.903	billion m3
The minimum release	1200	1200	m3/s
The minimum output	1000	1000	MW
Installed capacity	12,600	6000	MW
Efficiency coefficient	8.8	8.8	-

.

4.2. Data Description

The historical inflow of XLD hydropower reservoir is used as the input of the established optimal operation model of XLD and XJB hydropower reservoir system. Since the first unit of XLD hydropower station was put into operation in July, 2013, the historical inflow series are from 2014 to 2020. The optimization algorithm combined DP, POA, and DPSA is used to solve the established optimal operation model. Then the optimal operation process data of XLD and XJB hydropower reservoir system, including water level and outflow, are used as the input of optimal operation rule derivation model MVQIHGSO-SVM.

Specifically, the obtained daily operation data of seven years (from 2014 to 2020) of XLD and XJB hydropower stations were obtained with a total of 2557 observed sample data. Five years of the observed sample data from 2014 to 2018 were used as the training set of SVM model and the remaining data are used as validation set. Specifically, in order to improve the fitting accuracy, the end water level during a period was selected as decision variable. At the same time, in the implicit stochastic optimal operation of hydropower reservoirs, the initial water level and inflow of the period are often used as the input variables of the operation rules [9,56]. In addition, in order to consider the hydraulic connection between reservoirs, the initial water level of the same period of the adjacent reservoirs is also used as an input variables.

4.3. Results and Discussion

In order to verify the performance of the proposed MVQIHGSO-SVM model, other hybrid models including HGSO-SVM, PSO-SVM, SCA-SVM, and Grid-SVM models were compared in this study. The different validation indices of all the models are calculated according to Equations (33)–(36), which are given in

Table 8. Validation indices value of different operation rule derivation models on validation set.

Reservoirs	Models	R2	RMSE	MAE	MAPE
XLD	MVQIHGSO-SVM	0.998	0.340	0.126	0.021%
	HGSO-SVM	0.997	0.411	0.151	0.025%
	SCA-SVM	0.997	0.405	0.149	0.025%
	PSO-SVM	0.998	0.405	0.150	0.025%
	Grid-SVM	0.998	0.357	0.151	0.025%
XJB	MVQIHGSO-SVM	0.998	0.164	0.075	0.019%
	HGSO-SVM	0.997	0.192	0.098	0.026%
	SCA-SVM	0.997	0.189	0.092	0.025%
	PSO-SVM	0.996	0.192	0.096	0.025%
	Grid-SVM	0.997	0.187	0.086	0.023%

. Table 8 shows that for the same input data, the proposed MVQIHGSO-SVM model has more accurate results than the other models. Obviously, the equipped HGSO with multiple strategies cause a reduction in R², RMSE, MAE, and MAPE values in the prediction of end water level during a period. The best model according to the highest values of R², and the lowest values of RMSE, MAE and MAPE are the proposed MVQIHGSO-SVM as the best optimum model by a value of R² = 0.998, RMSE = 0.340, MAE = 0.126, and MAPE = 0.021% for XLD and R² = 0.998, RMSE = 0.164, MAE = 0.075, and MAPE = 0.019% for XJB. This means that the proposed MVQIHGSO method can find the best hyperparameter combination to improve the ability of SVM with respect to operation rule derivation. The remaining models, from best to worst according to RMSE value, are Grid-SVM, SCA-SVM, PSO-SVM, and HGSO-SVM for XLD and Grid-SVM, SCA-SVM, PSO-SVM, and HGSO-SVM for XJB. The worst model for XLD and XJB of all hybrid models are HGSO-SVM. The performance of the proposed model is improved by 17.27% and 14.58% compared with HGSO in term of RMSE value for XLD and XJB, respectively. The good performance of the proposed model in predicting data can be attributed to the following two reasons: (1) SVM shows unique advantages in solving small sample and non-linear problems by using kernel functions; (2) the multiple strategies including MVO algorithm and QI strategy improve the optimization ability of the basic HGSO.

In order to clearly show the effect of different methods to predict the end water level during a period, scatter plots of all hybrid models in the validation set are presented in Figure 4. Figure 4 shows that there is a significant quantitative correlation trend between observed and predicted water level process because the plotted points lie closely on a 45 degree straight line through the origin. However, it can be observed that performance of different hybrid models is slightly different and the proposed MVQIHGSO-SVM has a better performance than other hybrid models. For example, Figure 4 shows that the prediction error is small when XLD is at high water levels using MVQIHGSO-SVM, conversely, it is large using HGSO-SVM, SCA-SVM, and PSO-SVM. For XJB, the water level in the higher water level interval and the lower water level interval can be accurately predicted using MVQIHGSO-SVM model.

In order to further explore the difference in prediction accuracy of different hybrid models, the comparison of observed and predicted water level process of XLD and XJB are presented in Figure 5. As shown in Figure 5a for XLD, in the transitional stage from low water level to high water level, the prediction results of each hybrid model are quite different. The predicted water level process based on MVQIHGSO is closer to the observed values, followed by Grid-SVM, SCA-SVM, PSO-SVM, and HGSO-SVM. Here, combined with the inflow process of XLD from 1 January 2020 to 29 February 2020 in Figure 6, it can be observed that the inflow of XLD is relatively large after 23 January 2019. Therefore, the reservoir vacates the capacity in advance, resulting in the decline of the water level from 15 January 2019. From this point of view, the proposed MVQIHGSO-SVM can more accurately store and release floods in order to deal with future large and deterministic inflow processes. As shown in Figure 5b for XJB, when it is operating at a high water level, only the predicted water level process by MVQIHGSO-SVM meets the maximum water level constraint, namely the normal water level 380 m. This shows that the proposed MVQIHGSO-SVM model can more accurately find the operation rules of hydropower reservoirs considering the complex and non-linear constraints and the production experience of operators than other hybrid models.

Finally, the predicted water level processes based on all hybrid models are used as the input data of conventional operation model of XLD and XJB hydropower reservoirs. In addition, the optimal operation results of XLD-XJB are also calculated. These results are given in

Table 9. Total hydropower generation based on observed and predicted data as well as optimal operation.

Reservoirs	Hydropower Generation (TWh)
	Observed	MVQIHGSO-SVM	HGSO-SVM	SCA-SVM	PSO-SVM	Grid-SVM	Optimization
XLD	129.82	131.30	131.03	131.03	131.03	131.12	131.41
XJB	64.03	64.78	64.56	64.56	64.56	64.61	65.00
Total	193.85	196.08	195.59	195.58	195.58	195.73	196.41

. Table 9 shows that the highest total hydropower generation is obtained by optimal operation model, reaching 196.41 TWh. MVQIHGSO-SVM ranks second, followed by Gird-SVM, SCA-SVM, HGSO-SVM, PSO-SVM, and conventional operation. The hydropower generation obtained by conventional operation is 193.85 TWh, which is 2.56 TWh less than the optimal operation result. The hydropower generation calculated by the proposed MVQIHGSO-SVM model is closer to the optimal operation result, and the hydropower generation increased the most compared to conventional scheduling, reaching 2.23 TWh, increasing by 1.15%. In summary, the operation rule of hydropower reservoirs derived by the proposed MVQIHGSO-SVM model outperform other hybrid models, which is able to obtain total hydropower generation closer to the optimal dispatch model.

5. Conclusions

In this study, a novel operation rule derivation method combining improved HGSO and SVM is proposed. Multiple strategies including MVO operator and QI method are used to cope with the drawbacks faced by the original HGSO. Then the improved HGSO called MVQIHGSO is applied to optimize the hyperparameters of SVM model to derive the operation rules of hydropower reservoirs, forming MVQIHGSO-SVM model. The main contributions of this paper can be concluded below.

(1)

Multiple strategies are equipped into HGSO to improve its performance in exploration and exploitation. The multi-verse optimizer (MVO) is used to enhance the exploration capability of basic HGSO and help the inferior agent to escape from local optimal. Quadratic interpolation (QI) is used to improve the exploitation ability of HGSO. Finally, the exploration and exploitation are balanced by integrating the multiple strategies.

(2)

MVQIHGSO with multiple strategies is benchmarked by 23 classical benchmark functions. The results demonstrates that MVQIHGSO outperforms most of the well-known metaheuristic algorithms and has a superior efficacy compared to the competitors based on the convergence accuracy and speed.

(3)

MVQIHGSO-SVM model is used to derive operating rules of hydropower reservoirs. The XLD and XJB in the upper Yangtze River are selected as a case study. The results indicate that the proposed MVQIHGSO-SVM model can accurately obtain the joint operation rules of hydropower reservoirs. The total hydropower generation calculated by the proposed hybrid model is closer to the optimal operation result, and the hydropower generation increased the most compared to conventional scheduling, reaching 22.25 × 10⁸ kWh, increasing by 1.15%.

In the future, the improved HGSO algorithm can be combined with other machine learning method for deriving operation rules of huge hydropower reservoirs.