4.1.1. Objective Functions

With the advantages of low pollutant emission and high operational efficiency, hydropower now plays an increasingly important role in the energy system throughout the world [68–70]. As a result, the cascade hydropower reservoirs are often responsible for a variety of operational requirements. Considering the practical requirements of generation enterprises and electrical power systems in China, the generation benefit and peak operations were chosen as the focus of this paper.

(1) Power Generation

For almost all hydropower enterprises, the economic benefit is an important indicator of performance check under the market environment [71–73]. Generally, the same amount of available water resource can produce numerous scheduling schemes with different generation benefits. Hence, the first objective was chosen to find the best scheme that could maximize the total power generation of all hydropower reservoirs, which could be expressed as below:

$$E = \max \sum\_{t=1}^{T} \sum\_{n=1}^{N} P\_{n,t} \times \Delta\_t \tag{12}$$

where *E* is the power generation of all hydropower reservoirs. *T* is the number of periods. *N* is the number of hydropower reservoirs. *Pn*,*<sup>t</sup>* is the output of the *n*th hydroplant at the *t*th period. Δ*<sup>t</sup>* is the total hours of the *t*th period.

(2) Peak Operation

In modern power systems, the hydropower system is often asked to respond to the peak loads of power grid to produce a relatively smooth load series for other kinds of energy systems (like thermal, solar, and wind) [74–76]. In this way, the total operational cost of the electrical power system in the long run can be sharply reduced [77–79]. Hence, the second objective function was chosen to minimize the mean square deviation of the residual load curve obtained by subtracting all power outputs of hydropower reservoirs from the original energy load curve, which could be expressed as:

$$F = \min \sqrt{\frac{1}{2} \sum\_{t=1}^{T} \left[ Laad\_t - \sum\_{n=1}^{N} P\_{n,t} \right]^2} \tag{13}$$

where *Loadt* is the load demand of the power system at the *t*th period.

#### 4.1.2. Physical Constraints

(1) Water Balance Equation:

$$V\_{n,t+1} = V\_{n,t} + 3600\Delta\_t \times \left[ I\_{n,t} + \sum\_{i=1}^{\text{Num}} \left( Q\_{i,t} + S\_{i,t} \right) - Q\_{n,t} - S\_{n,t} \right] \tag{14}$$

where *Vn*,*t*, *In*,*t*, *Qn*,*t*, *Sn*,*<sup>t</sup>* are the storage volume, local inflow, turbine discharge and spillage of the *n*th hydroplant at the *t*th period, respectively. *Num* is the number of upstream hydroplants directly connected to the *n*th hydroplant.

(2) Reservoir Storage Constraint:

$$
\underline{V}\_{n,t} \le V\_{n,t} \le \overline{V}\_{n,t} \tag{15}
$$

where *Vn*,*<sup>t</sup>* and *Vn*,*<sup>t</sup>* denote the maximum and minimum storage volumes of the *n*th hydroplant at the *t*th period, respectively.

(3) Turbine Discharge Constraint:

$$
\underline{Q}\_{n,t} \le Q\_{n,t} \le \overline{Q}\_{n,t} \tag{16}
$$

where *Qn*,*<sup>t</sup>* and *Qn*,*<sup>t</sup>* denote the maximum and minimum turbine discharges of the *<sup>n</sup>*th hydroplant at the *t*th period, respectively.

(4) Total Discharge Constraint:

$$
\underline{q}\_{n,t} \le q\_{n,t} \le \overline{q}\_{n,t} \tag{17}
$$

where *qn*,*<sup>t</sup>* and *q n*,*t* denote the maximum and minimum total discharges of the *n*th hydroplant at the *t*th period, respectively.

(5) Power Output Constraint:

$$
\underline{P}\_{n,t} \le P\_{n,t} \le \overline{P}\_{n,t} \tag{18}
$$

where *Pn*,*<sup>t</sup>* and *Pn*,*<sup>t</sup>* denote the maximum and minimum power outputs of the *n*th hydroplant at the *t*th period, respectively.

(6) Initial and Final Water Levels Constraint:

$$\begin{cases} Z\_{n,0} = Z\_n^{\text{bcg}} \\ Z\_{n,T} = Z\_n^{\text{end}} \end{cases} \tag{19}$$

where *Zbeg <sup>n</sup>* and *Zend <sup>n</sup>* represent the initial and final water level of the *n*th hydroplant, respectively.

### *4.2. Technique for Order Preference by Similarity to an Ideal Solution (TOPSIS)*

As mentioned above, the multi-objective operation of cascade hydropower reservoirs is a typical multiple-criteria decision-making problem. Thus, it was necessary to find some tools to transform the original multi-objective problem into the single-objective optimization problem that can be addressed by EGSA. After a comprehensive comparison, TOPSIS was introduced to achieve this goal. In TOPSIS, the best scheme had the smallest distance from the positive ideal scheme and the largest distance from the negative ideal scheme, respectively. Then, the procedure of TOPSIS was give as below:

Step 1: Create an initial decision evaluation matrix with *m* schemes and *J* attributes. For the target problem, each scheme is composed of two attributes—the power generation *E* and the peak operation *F*. Then, the decision-making matrix *A* can be expressed as follows:

$$A = \begin{pmatrix} a\_{i,j} \end{pmatrix}\_{m \times f} = \begin{bmatrix} a\_{11} & a\_{12} & \cdots & a\_{1f} \\ a\_{21} & a\_{22} & \cdots & a\_{2f} \\ \vdots & \vdots & & a\_{i,f} & \vdots \\ a\_{m1} & a\_{m2} & \cdots & a\_{mf} \end{bmatrix} \tag{20}$$

where *ai*,*<sup>j</sup>* is the original value of the *j*th attribute in the *i*th scheme.

Step 2: Obtain the modified decision evaluation matrix to reduce the potential impact of various attribute properties (like measurement unit and tendency); given below:

$$B = \begin{pmatrix} b\_{i,j} \end{pmatrix}\_{m \times f} = \begin{bmatrix} b\_{11} & b\_{12} & \cdots & b\_{1f} \\ b\_{21} & b\_{22} & \cdots & b\_{2f} \\ \vdots & \vdots & b\_{i,j} & \vdots \\ b\_{m1} & b\_{m2} & \cdots & b\_{mf} \end{bmatrix} \tag{21}$$

$$b\_{ij} = \begin{cases} \frac{a\_j^{\max} - a\_{ij}}{a\_j^{\max} - a\_{\min}^{\min}} \text{ minimization attribute} \\\frac{a\_{ij} - a\_{\min}^{\min}}{a\_j^{\max} - a\_{\min}^{\min}} \text{ maximization attribute} \end{cases}, i = 1, 2, \cdots, m; j = 1, 2, \cdots, l \tag{22}$$

where *a*max *<sup>j</sup>* and *<sup>a</sup>*min *<sup>j</sup>* are the maximum and minimum values of the *j*th attribute. *Water* **2019**, *11*, 2040

Step 3: Obtain the normalized decision-making matrix to reduce the possible numerical problem caused by the feature differences (like units or magnitude), which could be expressed as below:

$$R = \begin{pmatrix} r\_{i,j} \end{pmatrix}\_{m \times f} = \begin{bmatrix} r\_{11} & r\_{12} & \cdots & r\_{1f} \\ r\_{21} & r\_{22} & \cdots & r\_{2f} \\ \vdots & \vdots & r\_{i,j} & \vdots \\ r\_{m1} & r\_{m2} & \cdots & r\_{mf} \end{bmatrix} \tag{23}$$

$$r\_{ij} = \frac{b\_{ij}}{\sqrt{\sum\_{i=1}^{m} b\_{ij}^2}}, i = 1, 2, \dots, m; \ j = 1, 2, \dots, J \tag{24}$$

where *rij* is the normalized value of the *j*th attribute in the *i*th scheme.

Step 4: Set the weight vector *<sup>w</sup>* = *w*1, ··· , *wj*, ··· , *wJ* , where *wj* is the weight of the *j*th attribute and there is *J j*=1 *wj*= 1. Then, calculate the weighted normalized decision-making matrix *V* by:

$$V = \left(v\_{i,j}\right)\_{m \times J} = \left(w\_j \times r\_{i,j}\right)\_{m \times J} \tag{25}$$

Step 5: Determine the positive ideal scheme *C*<sup>+</sup> and negative ideal scheme *C*<sup>−</sup> by:

$$\begin{cases} \mathbb{C}^{+} = \left[ c\_1^+, c\_2^+, \dots, c\_{\lceil \frac{l}{l} \rceil}^+ \right] = \max \{ v\_{i,j} | i = 1, 2, \dots, m \} \\ \mathbb{C}^{-} = \left[ c\_1^-, c\_2^-, \dots, c\_{\lceil \frac{l}{l} \rceil}^- \right] = \min \{ v\_{i,j} | i = 1, 2, \dots, m \} \end{cases} \quad j = 1, 2, \dots, l \tag{26}$$

Step 6: Calculate the Euclidean distances between the decision scheme and the positive (negative) ideal scheme, which could be expressed as below:

$$\begin{cases} D\_i^+ = \sqrt{\sum\_{j=1}^l \left( v\_{ij} - c\_j^+ \right)^2} \\\ D\_i^- = \sqrt{\sum\_{j=1}^l \left( v\_{ij} - c\_j^- \right)^2} \end{cases}, i = 1, 2, \dots, m \tag{27}$$

where *D*<sup>+</sup> *<sup>i</sup>* or *D*<sup>−</sup> *<sup>i</sup>* is the distance between the *i*th scheme and the position (or negative) ideal scheme. Step 7: Calculate the closeness of each decision plan to the negative ideal scheme by

$$f\_i = \frac{D\_i^-}{D\_i^- + D\_i^+}, i = 1, 2, \dots m \tag{28}$$

Step 8: Use the obtained closeness values to sort all the schemes, and the scheme possessing the maximum closeness will be treated as the optimal decision scheme.

$$f\_{\text{best}} = \max\{f\_1, f\_2, \dots, f\_m\} \tag{29}$$

#### *4.3. Details of EGSA for Multi-Objective Operation of Cascade Hydropower Reservoirs*

#### 4.3.1. Individual Structure and Swarm Initialization

To enhance the execution efficiency of reservoir operation, the total discharge of the hydropower reservoir is selected as the decision variable for encoding, and then the *i*th agent at the *k*th iteration *Xi*(*k*) can be expressed as below:

$$X\_i(k) = \begin{bmatrix} q\_{1,1} & q\_{1,2} & \cdots & q\_{1,T} \\ q\_{2,1} & q\_{2,2} & \cdots & q\_{2,T} \\ \vdots & \vdots & q\_{n,t} & \vdots \\ q\_{N,1} & q\_{N,2} & \cdots & q\_{N,T} \end{bmatrix} \tag{30}$$

During the initialization phase, the elements of the solution *Xi*(*k*) are randomly generated in the preset of the total discharge scopes of all hydroplants, which could be expressed as follows:

$$q\_{n,t} = \underline{q}\_{n,t} + random \times \left(\overline{q}\_{n,t} - \underline{q}\_{n,t}\right) \tag{31}$$

where *random* is the random number uniformly distributed in the range of [0,1].

#### 4.3.2. Heuristic Constraint Handling Method

During the evolutionary process, it is possible that some agents violate the equality or inequality constraints imposed on the hydropower system. In order to effectively modify the infeasible agents, this section presents a heuristic constraint handling method that tries to equally allocate the violated final storage volume to the turbine discharges of all adjusting periods. The detailed procedures for the solution *Xi*(*k*) are given as below:

Step 1: Set *n* = 1 and determine the necessary parameters for the constraint processing.

Step 2: Set the counter *L* = 1, and then use the water balance equation to calculate the storage of the *n*th hydroplant at the *t*th period by

$$V\_{n,t} = \max\left\{ \underline{V}\_{n,t}, \min\left\{ V\_{n,t-1} + 3600 \Delta\_{t-1} \times \left[ I\_{n,t-1} + \sum\_{i=1}^{\text{Num}} \left( Q\_{i,t-1} + S\_{i,t-1} \right) - Q\_{n,t-1} - S\_{n,t-1} \right] \overline{V}\_{n,t} \right\} \right\} \tag{32}$$

Step 3: Set *L* = *L* + 1, and then calculate the violation value η*<sup>n</sup>* of the final storage by Equation (33). Then, if η*n* is smaller than the accuracy <sup>ε</sup> or *<sup>L</sup>* is larger than the maximum iteration, go to Step 5; otherwise, turn to Step 4.

$$
\eta\_n = f\_n(Z\_{n,T}) - f\_n(Z\_n^{\text{cnd}}) \tag{33}
$$

where *fn*(·) is the nonlinear stage-storage curve of the *n*th hydropower reservoir.

Step 4: Use Equation (34) to obtain the modified total discharges of the target reservoir, and then calculate the corresponding reservoir storage volume by Equation (32) before turning to Step 3.

$$q\_{n,t} = \max\left\{\underline{q}\_{n,t}, \min\left\{q\_{n,t} + \frac{\eta\_n}{T \times 3600 \Delta\_t}, \overline{q}\_{n,t}\right\}\right\} \tag{34}$$

Step 5: Set *n* = *n*+1. If *n* ≤ *N*, turn to Step 2 for the new hydroplant; otherwise, stop the iteration.

#### 4.3.3. Calculation of the Modified Objective Values

Generally, the agents modified by the above heuristic constraint handling procedures would satisfy most of the physical constraints. Nevertheless, the agents might still violate some constraints due to reasons like unreasonable operation trajectory and limited adjusted times. To eliminate the negative effect of the infeasible agents, the value of the constraint violation involved in the solution *Xi*(*k*) (*Viol*[*Xi*(*k*)] for short) was obtained by Equation (35), and then merged into two original objective values (*E* and *F*) by Equation (36) to obtain the modified objective values (*E*' and *F*' ).

$$\operatorname{Void}[\mathbf{X}\_i(k)] = \sum\_{l\_1=1}^{L\_\mathcal{E}} c\_{l\_1} \cdot \max\left\{ \operatorname{g}\_{l\_1}(\mathbf{X}\_i(k)), 0 \right\} + \sum\_{l\_2=1}^{L\_\mathcal{E}} c\_{l\_2} \cdot \left[ \operatorname{e}\_{l\_2}(\mathbf{X}\_i(k)) \right]^2 \tag{35}$$

$$\begin{cases} E'[\mathbf{X}\_i(k)] = \sum\_{t=1}^T \sum\_{n=1}^N P\_{n,t} \times \Delta t - Viol[\mathbf{X}\_i(k)] \\\ F'[\mathbf{X}\_i(k)] = \sqrt{\frac{1}{2} \sum\_{t=1}^T \left[Load\_t - \sum\_{n=1}^N P\_{n,t}\right]^2} + Viol[\mathbf{X}\_i(k)] \end{cases} \tag{36}$$

where *gl*<sup>1</sup> (·) ≤ 0 is the *l*1th inequality constraint; *el*<sup>2</sup> (·) ≤ 0 is the *l*2th equality constraint; *cl*<sup>1</sup> and *cl*<sup>2</sup> are the penalty coefficients for the relevant inequality or equality constraint. *Lg* and *Le* denote the number of inequality and equality constraints, respectively.

4.3.4. Execution Procedures of the EGSA Method for the Target Problem

The flowchart of the developed EGSA method for solving multi-objective operation of cascade hydropower reservoirs is given in Figure 4.

**Figure 4.** The flowchart of enhanced gravitational search algorithm (EGSA) for multi-objective operation of cascade hydropower reservoirs.

#### **5. Case Studies**

#### *5.1. Engineering Background*

In this section, five reservoirs of the Wu hydropower system in Southwest China were chosen to test the performance of the proposed method, including Hongjiadu (HJD), Dongfeng (DF), Suofengying (SFY), Wujiangdu (WJD), and Goupitan (GPT). Since being put into operation, the Wu hydropower system has played a vitally important role in the environment protection and economic development of the Guizhou province. Figure 5 shows the topological structure of the selected hydropower reservoirs, while the typical load demand curves of the four seasons are shown in Figure 6. From Figures 5 and 6, it could be observed that in the Wu hydropower system, there were tight hydraulic and electrical connections from upstream to downstream reservoirs, while four load curves with different features further increased the optimization difficulty.

**Figure 5.** The topological structure of the studied hydropower reservoirs.

**Figure 6.** Typical load curves in different season.

#### *5.2. Case Study 1: Power Generation of Cascade Hydropower Reservoirs*

#### 5.2.1. Robustness Testing of Different Evolutionary Algorithms

In this section, to verify the robustness of the developed method, four famous methods (DE, PSO, SCA, and GSA) were introduced for comparative analysis. Table 5 shows the detailed statistical results of the 5 methods in 20 independent runs for 4 cases, where the swarm size and maximum iteration were set as 50 and 500, respectively. It can be clearly seen that the solutions obtained by the EGSA method were more stable than the other methods. For instance, in Case 1, the EGSA method could make about 96.4%, 96.5%, 97.8%, and 97.7% reductions in the standard deviations of the DE, PSO, SCA, and GSA. Hence, the feasibility of the EGSA method was fully demonstrated in this case.

**Table 5.** Statistical results of the 5 methods in 20 independent runs for 4 cases of power generation (104 kW·h).


Figure 7 draws the box and whisker test of the best-so-far solutions obtained by 5 algorithms in 4 cases. It can be observed that the performances of four control methods in the power generation of hydropower system are relatively stable but still inferior to the EGSA method possessing the smallest variation ranges in the final solutions. Therefore, the proposed method proves to be an effective tool to deal with the power generation of a hydropower system.

**Figure 7.** Box and whisker test of different algorithms in different cases for power generation.

#### 5.2.2. Comparison of the Optimal Results Obtained by Different Evolutionary Algorithms

Table 6 shows the detailed results in the best solutions obtained by different methods. It can be seen that in all cases, EGSA can always find the best results among all algorithms. For instance, compared with DE, PSO, SCA, and GSA, the EGSA method could make about 8.77 <sup>×</sup> 104 kW·h, 3.2 <sup>×</sup> <sup>10</sup><sup>4</sup> kW·h, 8.61 <sup>×</sup> 10<sup>4</sup> kW·h, and 8.5 <sup>×</sup> 104 kW·h improvements in power generation in Case 1, respectively. Additionally, among all hydroplants, the GPT hydroplant with a huge installed capacity produced the largest proportion of power generation in all solutions, which was in line with the actual operation situation of the hydropower system [80–83]. Thus, the above analysis demonstrated the effectiveness and rationality of the scheduling process obtained by the proposed method.


**Table 6.** Detailed results in the best solutions by different methods (104 kW·h).


**Table 6.** *Cont*.

5.2.3. Convergence Analysis of Different Evolutionary Algorithms

Figure 8 shows the convergence trajectories of 5 algorithms in different cases. It could be observed that in four cases, SCA failed into a local optima at the early stage, and then started to improve the quality of solutions as the iteration number reached 300, but still failed to find out the satisfying solutions at the end; due to the loss of individual diversity [84], three other methods (PSO, GSA, and DE) outperformed the SCA method but were still inferior to the proposed method from the beginning to end. Additionally, the proposed method could quickly converge to the scheduling schemes that were close to the optimal solution at the early stage, and then change slightly with an increasing number of iterations. Thus, it could be concluded that three modified strategies could effectively enhance the performance of the standard GSA method.

**Figure 8.** Convergence trajectories of 5 methods in different cases.

*5.3. Case Study 2: Peak Operation of Cascade Hydropower Reservoirs*

5.3.1. Robustness Testing of Different Evolutionary Algorithms

In this section, the proposed method was applied to deal with the peak operation of a hydropower system. Table 7 shows the statistical results of the 5 methods in 20 independent runs for 4 cases. It was observed that all of the methods could achieve the goal of reducing the peak loads of a power system due to the reductions in the statistical indices; while the proposed method had the performances among all methods. For instance, the average improvements in the range of the objective values of EGSA was about 99.7%, compared to the three other methods in the 4 cases. Hence, the feasibility of the EGSA method in addressing the peak operation of the hydropower system was proved in this case.


**Table 7.** Statistical results of the 5 methods in 20 independent runs for 4 cases of peak operation.

Figure 9 shows the box and whisker test of 5 algorithms in different cases. It can be seen that the distribution ranges of the solutions obtained by the proposed method were much smaller than the other methods, demonstrating the effectiveness of the constraint handling method and the modified strategies. Thus, it can be concluded that the EGSA method could produce stable scheduling schemes for the complex hydropower system peak operation problem.

**Figure 9.** Box and whisker test of 5 algorithms in different cases for peak operation.

#### 5.3.2. Comparison of the Optimal Results Obtained by Different Evolutionary Algorithms

To clearly show the feasibility of the EGSA method, the detailed results in the best solutions as obtained by different methods are given in Table 8. It was found that the EGSA method could achieve satisfactory results. For instance, the proposed method could bring about a 30.9% reduction in the peak values in Summer, which was obviously better than 21.0% of DE and 24.9% of PSO, 23.9% of SCA and 19.4% of GSA. Thus, this simulation clearly proved the superiority of the proposed method in responding to the peak operation requirement of the power system.


Note: Reduction = original–optimization of the target method.

5.3.3. Rationality Analysis of the Best Results Obtained by the Different Evolutionary Algorithms

Figure 10 shows the detailed output process obtained by the different methods in four cases. It was found that there were often two peak periods (morning and evening) in the original load curves, and obvious differences in the load features (like peak times or numbers) of the four load curves, which would increase the optimization difficulty of the peak operation. The two methods (GSA and EGSA)

exhibited excellent responses to the load changes by collecting the hydropower generation at the peak periods and storing the energy at valley periods. The EGSA method was superior to the GSA method due to its smoother residual load curves; besides, the outputs of all the hydroplants varied in the feasible range between the minimum output limit and the installed capacity. Thus, the EGSA method could produce feasible solutions for the peak operation of a hydropower system.

**Figure 10.** Detailed output process obtained by the different methods in the four cases.

5.4.1. Comparative Analysis of the Optimal Results Obtained by the Different Methods with 100 Weights

The focus of this section is on the demonstration of the feasibility of the proposed method in addressing the multi-objective optimization problems of a hydropower system. Figure 11 draws the distributions of the objective functions obtained by different methods, where each method was executed in 100 different weight combinations. Specifically, the weight *w*<sup>1</sup> for power generation was increased from 0 to 1.0 at the same interval of 0.01, while the weight *w*<sup>2</sup> for peak operation was set as 1.0 – *w*1. From Figure 11, it was observed that the total generation was gradually decreasing with the increasing objective value of peak operation, which implied that there was an obvious conflict between power generation and peak operation. Additionally, the solutions of PSO and GSA were dominated by that of the EGSA method, which meant that the EGSA method had a higher probability to obtain the Pareto optimal solutions than the PSO and GSA method. Thus, the proposed method can generate the near-optimal Pareto solutions to approximate the relationship between power generation and peak shaving in practice.

*<sup>5.4.</sup> Case Study 3: Mutli-Objective Operation of Cascade Hydropower Reservoirs*

**Figure 11.** Distributions of the objective functions obtained by different methods with 100 weights.

5.4.2. Rationality Analysis of the Best Results Obtained by the Different Evolutionary Algorithms

Figure 12 shows the detailed results of three typical schemes obtained by the EGSA method in Spring, where Scheme 1 showed the maximum power generation, Scheme 3 had the best peak operation performance, while Scheme 2 could achieve a balance between power generation and peak operation. It could be seen that for all the hydropower reservoirs, the water levels varied in the preset range between a dead water level and a normal water level while the power outputs were smaller than the installed capacity, demonstrating the effectiveness of the heuristic constraint handling method in guaranteeing the feasibility of the solution. Meanwhile, there were obvious differences in the peak operation and power generations of the three typical schemes, which indicated that it was necessary to choose a suitable scheme based on the actual requirement. Thus, this case again proved the practicability of the EGSA method in solving the multi-objective operation problems.

**Figure 12.** *Cont*.

**Figure 12.** Detailed scheduling results of three typical schemes obtained by EGSA in Spring.

#### **6. Conclusions**

In this research, an enhanced gravitational search algorithm coupled with the Technique for Order Preference by Similarity to an Ideal Solution (TOPSIS) was developed to deal with the complex multi-objective operation problem of cascade hydropower reservoirs balancing the power generation and peak operation benefits. The proposed method and several famous evolutionary methods were used to solve the 12 famous benchmark functions and the optimal operation of the Wu hydropower system in China. The critical findings obtained from the experimental results are given as below:


**Author Contributions:** All authors contributed extensively to the work presented in this paper. Z.-k.F. and W.-j.N. contributed to the finalized the manuscripts. S.L. contributed to the modeling and programming. Z.-q.J. contributed to the literature review. M.-s.M. and B.L. contributed to the data analyses.

**Funding:** This research was supported by the Natural Science Foundation of Hubei Province (2018CFB573), National Natural Science Foundation of China (51709119 and 51909133), and the Fundamental Research Funds for the Central Universities (HUST: 2017KFYXJJ193).

**Acknowledgments:** The writers would like to express appreciation to both editors and reviewers for their valuable comments and suggestions.

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

#### **References**


© 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).
