*3.3. Normalization of Objective Sequence*

In this study, the three target values are normalized. To eliminate the difference of dimension of each objective and make them comparable, each object is normalized according to the following formula:

$$y\_{i\bar{j}} = \left[\mathbf{x}\_{i\bar{j}} - \min(\mathbf{x}\_{i\bar{j}})\right] / \left[\max(\mathbf{x}\_{i\bar{j}}) - \min(\mathbf{x}\_{i\bar{j}})\right] \tag{10}$$

where, *xij* and *yij* are target eigenvalue and the normalized value of the non-inferior solution *j*(*j* = 1, 2, ... , *n*) of the object *i*(*i* = 1, 2, ... , *m*) respectively; Δ*xij* is the difference between the target eigenvalue and the moderate value.

Among the three objectives of this study, the average irrigation water shortage (Obj-1) and ecological water shortage (Obj-2) are the minimizing objectives. Conversely, the multi-year average generation capacity (Obj-3) belongs to the maximizing objective. The optimization directions of objectives are different. In order to construct the joint sequence values of two targets, the optimization direction of normalized values needs to be consistent. The normalized value of Obj-1 and Obj-2 were subtracted by 1. A new normalized matrix is generated by replacing the original normalized value with the difference.

#### *3.4. Copulas Theory*

The copula was originally introduced by Sklar et al. [40] as a useful method to derive joint distributions with marginals. Its design allows it to deal with non-normal distributions. This method is also likely to be beneficial when the marginals can be definitely stipulated but the joint distribution is not straightforward to construct. The name 'copula' basically comes from the word 'couple' to emphasize a manner whereby a joint distribution function and the marginals can be combined [23]. The copula is a multidimensional joint probability distribution function with a uniform distribution in [0,1] interval. It can connect multiple random variables *F*(*x*1, *x*2, ... , *xn*) with their respective marginal distribution functions *F*1(*x*1), *F*2(*x*2), ... , *Fn*(*xn*) for the joint distribution function. Therefore, copula is often referred to as 'connection functions' or 'dependent functions.' Among them, the marginal distribution function describes the distribution of a single variable, while the copula function indicates the joint distribution of multi-dimensional composite variables [41].

The 'Sklar theorem' defines that: Assuming *F* is an *n*-dimensional distribution function whose marginals are (*F*1, *F*2, ... , *Fn*) of a random vector (*X*1, *X*2, ... , *Xn*), there exists an *n*-dimensional copula function *C*, which satisfies the following formula for arbitrary *x* ∈ *Rn*:

$$F(\mathbf{x}\_1, \mathbf{x}\_2, \dots, \mathbf{x}\_n) = \mathbb{C}[F\_1(\mathbf{x}\_1), F\_2(\mathbf{x}\_2), \dots, F\_n(\mathbf{x}\_n)] = \mathbb{C}(\mathbf{u}\_1, \mathbf{u}\_2, \dots, \mathbf{u}\_n) \tag{11}$$

If the marginal distributions *F*1, *F*2, ... , *Fn* are continuous, function *C* is uniquely determined. Conversely, if *C* is an *n*-dimensional copula function and *F*1, *F*2, ... , *Fn* are a set of univariate distribution functions, the function *F* defined by Equation (11) is an *n*-dimensional distribution function having marginal distribution as *F*1, *F*2, ... , *Fn*.

Copula functions could be divided into four categories as a whole: Archimedean Copula, Extreme Copula, Elliptical Copulas and other hybrid families of copulas. In this study, two-dimensional copula of Clayton Copula, Frank Copula and Gumbel Copula in the Archimedean Copula function family and Gaussian Copula and Student Copula in the elliptic population are selected.

Copula is able to construct multi-dimensional joint probability distribution function by arbitrary marginal distribution and correlation structure [42]. It is flexible in forms and does not require all variables to obey the same scatter. In recent years, copula is often used to analyze the joint recurrence period of multi-variables and the frequency of combined events. Copula function has

natural advantages in constructing joint distribution of two variables. In this study, copula is used to construct joint sequence values of two targets and applied to the study of multi-objective competition mechanism for the first time. Procedures of constructing joint sequence values of two targets are as follows:


The AIC (Akaike Information Criterion) method is proposed by Akaike in 1974 [46] from the perspective of information theory. It is a goodness-of-fit evaluation criterion based on the Kullback-Leibler information metric, which contains two factors: the deviation between the empirical points and theoretical copula functions as well as the error fluctuation caused by the number of parameters of copula functions. Taking the joint copula function of two variables as an example, the concrete formulas of the AIC information criterion method are as follows:

$$MSE = \frac{1}{n} \sum\_{i=1}^{n} \left( F\_n(x\_{i\cdot} y\_i) - \mathbb{C}(u\_{i\cdot}, v\_i) \right)^2 \tag{12}$$

$$AIC = n \ln(MSE) + 2m \tag{13}$$

where, *Fn*(*xi*, *yi*), *C*(*ui*, *vi*) are the empirical and theoretical frequencies of the joint copula function with two variables are expressed, respectively. *MSE* is the mean square error and *m* is the number of parameters of the model.

Among them, the smaller the value of AIC, the better the fitting of the copula function. The AIC information criterion is applicable to the comparison and optimization of copula functions with a different number of model parameters.

5. Computation of Joint Distribution Sequence Values. In order to find a new sequence that reflects the overall characteristics of the two variables, it needs to inverse the frequency sequence obtained above and then get the joint sequence value of the two variables. In this study, the inverse function, that is, NORMINV for the normal cumulative distribution function in the MATLAB software is used to derive the sequence values of the joint probability distribution of the copula function.

In this paper, with the Pareto set output from the multi-objective optimal operation model as the input, copula is employed to study the multi-objective competition mechanism. The main research process is as follows (Figure 2):

**Figure 2.** The main research framework.

#### **4. Results and Discussions**

## *4.1. Optimal Results of the ICGC-NSGA-II*

In the main stream of the Heihe River, Huangzangsi Reservoir, seven run-off hydropower stations of Baopinghe, Sandaowan, Erlongshan, Dagushan, Xiaogushan, Longshou-II and Longshou-I and Zhengyixia Reservoir at the end of the middle reaches constitute the structure of '2 reservoirs 7 stations' as the cascade reservoirs' hydropower station system. With the ICGC-NSGA-II algorithm for solving the optimal dispatching multi-objective model, the decision variables are the final water level of the time interval of Huangzangsi Reservoir and Zhengyixia Reservoir. The population size is 2000 and the number of iterations is 500 generations. The results of long-term optimal scheduling are sorted according to the ascending order of annual water shortage, as shown in Table 1.

According to Table 1, the minimum and maximum annual average water shortage of irrigation are 925 <sup>×</sup> 104 m3 and 6595 <sup>×</sup> 104 m3; the minimum and maximum of the ecological water shortage are 3925 <sup>×</sup> 104 m3 and 9058 <sup>×</sup> 10<sup>4</sup> m3; and the minimum and maximum of the multi-year average power generation are 25.54 <sup>×</sup> 108 kW·h and 26.99 <sup>×</sup> 10<sup>8</sup> kW·h, respectively. The non-inferior solution set is plotted in one coordinate system and the Pareto surface maps of irrigation, ecology and power generation are obtained, as shown in Figure 3.


**Table 1.** Multi-objective Pareto optimal set.

**Figure 3.** Multi-objective Pareto surface graph of Obj 1~3 in one coordinate system. Blue points indicates non-inferior solution.

Intuitively, the three-dimensional spatial distribution and the surface shape of the multi-objective non-inferior solution are shown in Figure 3. The non-inferior solutions appear in clusters in the Pareto surface graph (Figure 3) and are distributed over a wide range. The results show that there is little difference among non-inferior solutions in the same cluster and the optimal scheduling results converge well near the clusters to which the non-inferior solution belongs. The algorithm itself retains the diversity of non-inferior solutions. As the basic regulating rules, the Obj-1 and Obj-2 are 'the smaller the better' and the Obj-3 is the larger the better. Hence, theoretically, the optimal solution for the scheme is '925 <sup>×</sup> 104 m3, 3925 <sup>×</sup> 104 m3, 26.99 <sup>×</sup> 108 kW·h.' However, the three objectives cannot reach the optimal solution of the sample theory at the same time—an indication of mutual restrictions and influences among the three objectives. The Pareto surface of the scheme is a smooth space surface oriented to the vector ± (−1, −1, 1) direction (the main vertical line of the surface is consistent with the direction of a vector, which is expressed as the surface oriented to the vector direction). It can be inferred that the set of non-inferior solutions is composed of the non-inferior solution, which is weighted by the equivalent distance from the optimal solution of the sample theory. Vector ± (−1, −1, 1) is the superior direction of non-inferior solutions, consistent with the direction of the non-inferior solutions set towards the optimal solution of the sample theory, while vector ± (1, 1, −1) is the inferior direction.

It is difficult to exhibit the mutual restriction among different objectives; therefore, the three-dimensional stereogram was converted into two-dimensional plane map [47] and three objectives are respectively described in terms of bubble size. The bigger the bubble, the better the objective. Thus, the multi-objective bubble diagram is generated and shown in Figures 4–6.

**Figure 4.** Bubble graph of irrigation and ecology, in which the bubble size represents the multi-year average power generation (Obj-3).

**Figure 5.** Bubble graph of irrigation and power generation, in which the bubble size represents the multi-year average ecological water shortage (Obj-2).

**Figure 6.** Bubble graph of ecology and power generation, in which the bubble size represents the multi-year average irrigation water shortage (Obj-1).

The bubble diagram shows an approximate linear distribution with layers (Figures 4–6). Taking the bubble graph of irrigation and ecology (Figure 4) as an example, the majority distribution of non-inferior solutions spreads wider. The stratification phenomenon and the wide distribution are supportive of the convergence and diversity of the optimal dispatching results. The approximate linear distribution indicates that there exists a distinct growth-decline phenomenon and a direct water competition relationship between the irrigation water shortage and the ecological water shortage. It is notable that the closer to the intersection of the two coordinate axes, the smaller the bubble; while the farther away from the intersection, the larger the bubble. Within the same layer, the bubble size is equivalent. According to the basic regulating rules, it is evident that the closer the equivalent weighted distance from the intersection, the better the non-inferior solution and the larger the air bubbles. The distribution of bubble size illustrates the competitive relationship between power generation with irrigation and ecological water, that is, the smaller the water shortage of irrigation and ecology, the less power generation capacity. The contrary indicates more power generation.

#### *4.2. Analysis of Two-Objective Competition Mechanism*

As a way of exploring the competition mechanism and quantitatively describe the transformation law among objectives, the Pareto frontier is fitted to obtain the quantitative transformation formula between any two objectives. According to the previous analysis, the three-objective Pareto solution set is fully optimized and most of the solutions are high-quality feasible solutions close to the optimum one. On the premise of sufficient optimization, there is a macro-rule of 'one falls another rises.' By sorting the normalized values of one objective and drawing two-dimensional scatter plots of the normalized indices of the other two, several solutions are intercepted using the normalized values from small to large. It can be found that these solutions are exactly the Pareto frontier of the other two objectives. In other words, when the normalized value of one objective is small (that is, the solution is inferior), then the other two targets show the strongest regularity and optimum. Additionally, the number of interception solutions for sorting targets depends on the coverage of the Pareto frontier of the other two targets.

As shown in Figure 7, the two-dimensional scatter plots of Obj-1 and Obj-2 normalization indices, the non-inferior solutions are sorted in descending order according to the normalized value of Obj-3. After 700 solutions are intercepted, the Pareto frontier of Obj-1&2 is obtained as the orange points (Figure 7). A mutation point (0.7968, 0.4905) is found in the red dot (Figure 7).

**Figure 7.** The two-dimensional scatter plots of Obj-1 and Obj-2's normalized value. All the points indicate Pareto set, and among them, the orange points represent the Pareto frontier.

The scatter plots are drawn by selecting two sequences before and after the mutation points in the Pareto frontier (Figure 7). Then, the piecewise linear function relationship of the two series is obtained, as shown in Figure 7. The correlation coefficients of the two series are 0.981 and 0.942, respectively.

From the above formula, it can be seen that both the slopes of the two piecewise functions are less than zero and the normalized values of Obj-1 and Obj-2 are inversely correlated.

In front of the mutation point, the independent variable fluctuates within the range of [0,0.7968] and the normalized value of Obj-2 decreases by 0.7726 for each increase of Obj-1. Following the mutation point, the range changes to [0.7968,1] and the normalized value of Obj-2 decreases by 1.7388 for each increase of Obj-1. Converting to their own respective target values of the scheme, when the independent variable is between 925~2078 <sup>×</sup> <sup>10</sup><sup>4</sup> m3, the average annual ecological water shortage will increase by 9.0553 <sup>×</sup> 104 m<sup>3</sup> for each reduction of 1 <sup>×</sup> 104 m<sup>3</sup> of the average annual irrigation water shortage. Yet, when the independent variable is between 2078~6595 <sup>×</sup> 104 m3, the average annual ecological water shortage will increase by 0.6825 <sup>×</sup> 104 m3 for each reduction of 1 <sup>×</sup> 104 m3 of the average annual irrigation water shortage. Conclusively, in the interval where the irrigation water shortage is large, the increase of the ecological water shortage is slow while reducing the irrigation water shortage. In the interval where the irrigation water shortage is small, continuing to reduce the irrigation water shortage will lead to a significant increase in the ecological water shortage.

The two-dimensional scatter plot of the normalized indices of 'Obj-1 and Obj-3' and 'Obj-2 and Obj-3' is shown in Figures 8 and 9 in the way similar to Figure 7, with correlation coefficients of 0.984 and 0.933, respectively.

From Figure 8, it can be seen that the slope of the trend line is negative and the normalized values of Obj-1 and Obj-3 are inversely correlated. The normalized value of Obj-3 decreases by 1.3009 units for each additional unit of Obj-1 normalized value. An indication is that for each reduction of the average annual irrigation water shortage by 100 <sup>×</sup> 104 m3, the average annual power generation decreases by 0.03349 <sup>×</sup> 10<sup>8</sup> kW·h. The average annual irrigation water shortage is positively correlated with the average annual power generation.

From Figure 9, the slope is also negative and the normalized values of Obj-2 are inversely correlated with Obj-3. The normalized value of Obj-3 decreases by 1.0687 units for each additional unit of Obj-2 normalized value. The inference of the results is that: for each reduction of the average annual ecological water shortage by 100 <sup>×</sup> 104 m3, the average annual power generation decreases by 0.03443 <sup>×</sup> 108 kW·h. There is a positive correlation between the average annual irrigation water shortage and the average annual power generation.

**Figure 8.** The two-dimensional scatter plots of Obj-1 and Obj-3's normalized value. All the points indicate Pareto set, and among them, the orange points represent the Pareto frontier.

**Figure 9.** The two-dimensional scatter plots of Obj-2 and Obj-3's normalized value. All the points indicate Pareto set, and among them, the orange points represent the Pareto frontier.

#### *4.3. Analysis of Three-Objective Competition Mechanism*

Considering that the objectives of the Heihe River Basin are integrally and mutually restrictive, it is virtually impossible to reveal the law of the whole water resources system only by studying the relationship between the two objectives. This section synthesizes the sequence values of any two targets into a new sequence that represent the characteristics and information of two sub-sequences and analyses impact of the change of one target on the whole of the other targets in the Heihe River Basin. Through these analyses, the competition mechanism among the three objectives is explored.

Dependence between hydrological series is the premise of constructing the joint distribution using multivariate copula. The previous study indicates that the correlation coefficient between any two objectives is within 0.933~0.984. Such a high correlation supports the construction of the joint copula function of two variables.

The marginal distribution of each target is constructed with the non-parametric method of Gringorten and the empirical frequency estimates of each target sequence are obtained. The maximum likelihood method is used to get the parameter estimates of five joint copula functions between the two target sequence values and the corresponding cumulative distribution function values are calculated. The goodness of fit of five copula functions is evaluated with the AIC method which are best matching with existing hydrological sequences. The AIC evaluation indices of five combined copula functions under different objective combinations are obtained as shown in Table 2.


**Table 2.** Five copula functions' Akaike Information Criterion (AIC) evaluating value.

According to the AIC information criterion, the smaller the value of AIC evaluation index, the better the fit of representative copula function. As can be seen from Table 2 above, the minimum AIC evaluation index of Obj-1&2 is −16,437.86 with the Student Copula function. Similarly, the minimum in Obj-1&3 and Obj-2&3 are −10,671.15 (Gaussian Copula) and −11,292.73 (Frank Copula) respectively. As a result, the three two-variable-copula functions are selected as the optimal bivariate copula joint distribution for combining sequences.

To derive the sequence value of joint distribution, the joint sequence value of the optimal copula function under three combinations is obtained. The joint sequence value covers the characteristics and the information of the two objective sequences that represent the overall level. The scatter plot of the Obj-1 normalized value and the Frank Copula joint sequence value of Obj-2&3 is illustrated in Figure 10, as well as the 'Obj-2 + Obj-1&3(Gaussian Copula)' and 'Obj-3 + Obj-1&2(Student Copula)' showed in Figures 11 and 12, respectively.

Comparing Figures 10–12, Figure 12 appears the most discrete, which indicates that the regularity of the impact of Obj-3 on Obj-1&2 in the Heihe River Basin is the worst. It is known that power generation does not consume but rather utilizes surface water resources, so the competition between power generation and the other objectives is minimal. Moreover, the slope is the smallest, indicating that the change of power generation has the least impact on Obj-1&2.

**Figure 10.** The normalized value of Obj-1 and Obj-2&3's joint series scatter diagram.

**Figure 11.** The normalized value of Obj-2 and Obj-1&3's joint series scatter diagram.

**Figure 12.** The normalized value of Obj-3 and Obj-1&2's joint series scatter diagram.

The increase of Obj-1 normalization index will lead to the decrease of combined sequence value of Obj-2&3, which indicates that when Obj-1 tends to be more optimized, the Obj-2&3 will become worse as a whole. The slope of the Obj-1 normalized value in scatter plot ahead of 0.9 is milder than that after 0.9, which denotes that when irrigation water shortage is greater than 1492 <sup>×</sup> 104 m3, with the decrease of irrigation water shortage, the overall impact on ecological water and power generation is smaller and the cost of the optimal irrigation water shortage is lower. When Obj-1 is more satisfied, that is, when irrigation water shortage is less than 1492 <sup>×</sup> <sup>10</sup><sup>4</sup> <sup>m</sup>3, with the decrease of irrigation water shortage, the overall impact on ecological water and power generation is greater and the cost of the optimal irrigation water shortage is higher. Therefore, it is recommended that the average annual water shortage for irrigation should be about 1492 <sup>×</sup> <sup>10</sup><sup>4</sup> <sup>m</sup>3.

The plot of the Obj-1 normalization index is scattered before 0.6 and clustered after 0.6. It illustrates that when the irrigation water shortage is greater than 3193 <sup>×</sup> <sup>10</sup><sup>4</sup> <sup>m</sup>3, the change of irrigation water has little effect on the ecological water and power generation integrally; after 0.6, when the irrigation water shortage is less than 3193 <sup>×</sup> 104 <sup>m</sup>3, the irrigation water exhibits a strong influence on other objectives.

The increase of the Obj-2 normalization index will lead to the decrease of the combined sequence value of Obj-1&3, which indicates that when the Obj-2 is more optimized, the Obj-1&3 may become worse as a whole. The slope of the Obj-2 normalization index in the scatter plot before 0.8 is milder than that after 0.8, which denotes that, in the stage of lower ecological satisfaction, when ecological water shortage is greater than 4951 <sup>×</sup> <sup>10</sup><sup>4</sup> m3, with the decrease of ecological water shortage, the overall impact of irrigation water and power generation is smaller and the cost of optimizing the ecological water shortage is lower. Otherwise, when ecological water shortage is less than 4951 <sup>×</sup> 10<sup>4</sup> m3, with the decrease of ecological water shortage, the overall impact on irrigation water and power generation is greater and the cost of optimizing ecology water shortage is higher. It is recommended that the average annual water shortage in the ecological process should be around 4951 <sup>×</sup> 104 <sup>m</sup>3.

The increase of normalized index of Obj-3 will lead to a decrease of the combined sequence value of Obj-1&2, which indicates that when the Obj-3 tends to be more optimized, the Obj-1&2 will become worse as a whole. The scatter plots before 0.65 are discrete and after 0.65 are concentrated, showing that when the average annual power generation is less than 26.48 <sup>×</sup> 108 kW h, the change of power generation has little effect on the irrigation and ecological water as a whole; on the contrary, the impact is stronger.

#### **5. Conclusions**

This study takes the multi-objective joint optimal dispatch of cascade reservoirs in the Heihe River Basin as a study object. Based on the ICGC-NSGA-II algorithm to solve this model, the Pareto non-inferior solution set is obtained. The competition mechanism among two objectives is quantitatively described by statistical means and with the copula function constructing the joint sequence of two targets, the three objectives' competition mechanism is explored too. This study expects to provide a foundation for selective preference of the Pareto set and a new idea for multi-objective study. The main conclusions are summarized as follows:


To summarize, the copula function could combine the marginal distribution of any two sequences and construct a new joint sequence and all the information of the sub-sequence is contained, so there is no information distortion during the combination process. It is an effective tool for quantitatively studying the multi-objective competition mechanism. At present, the research on multi-objective competition mechanisms is in the preliminary stage and the methods adopted in this paper enrich this field.

**Author Contributions:** Conceptualization, software, data curation and writing—original draft preparation, M.Z.; methodology, L.W.; validation, H.W.; formal analysis, G.L.; resources, S.H.; writing—review and editing, S.L.; supervision, project administration and funding acquisition, Q.H.

**Funding:** This study was jointly funded by the National Key Research and Development Program of China (grant number 2017YFC0405900), the National Natural Science Foundation of China (grant number 51709221), the Planning Project of Science and Technology of Water Resources of Shaanxi (grant numbers 2015slkj-27 and 2017slkj-19), the China Scholarship Council (grant number 201608610170) and the Doctorate Innovation Funding of Xi'an University of Technology (grant number 310-252071712).

**Conflicts of Interest:** The authors declare no conflict of interest.
