Causal Inference of Optimal Control Water Level and Inﬂow in Reservoir Optimal Operation Using Fuzzy Cognitive Map

: Reservoir optimal operation (ROO) has always been a hot issue in the ﬁeld of water resources management. Analysis of the relationship of optimal control water level and inﬂow is conducive to understanding and solving ROO under deterministic inﬂow conditions. The current research uses a fuzzy cognitive map (FCM) as a tool to e ﬀ ectively model complex systems and then extracts systematic relationship diagrams from the dataset. A new fuzzy cognitive map with o ﬀ set (FCM-O) is proposed to overcome the causal inference error caused by non-linear mapping of the activation function in a traditional FCM. With the application of inferring the causal relationship between the optimal control water level and inﬂow of ROO for the Three Gorges Reservoir (TGR), the experimental results show that, compared with FCM in the min data error, FCM-O reduces 11.11% and 7.14% in the training and the testing, respectively. Also, the experimental results of FCM-O are more reasonable than those of FCM. Finally, the following conclusions about the causal inference of optimal control water level and inﬂow in ROO for TGR are drawn: (1) The optimal control water level in September, October and November needs to be raised as much as possible to raise the water head of power generation, which is mainly a ﬀ ected by the constraints of the maximum operating water level of the reservoir rather than inﬂow; (2) the optimal control water level in January, February and March is positively a ﬀ ected by the inﬂow of the adjacent months; (3) the optimal control water level in April is due to the approaching ﬂood season. In order to prevent water discarding, the water level is low and the optimum operation space is small. All of those shows that FCM-O is more competent than FCM in the causal relationship between optimal control water level and inﬂow in ROO.


Introduction
Reservoir operation optimization (ROO) facilitates flood control, agriculture irrigation, hydropower generating and shipping [1,2], which serves mankind by optimizing benefits through meeting societal demands [3]. ROO has always been a hot issue in the field of water resources management; many researchers have carried out multiple studies on reservoir optimal operation, such as: Power generation optimal operation [4][5][6][7], flood control optimal operation [8][9][10][11], multi-objective optimal operation [12][13][14][15][16][17][18], stochastic scheduling [19][20][21][22][23], and so on [24][25][26]. To solve ROO efficiently, traditional mathematical programming [27][28][29][30][31] and modern heuristic algorithms [32][33][34][35] are proposed. In addition, some studies on the characteristics of ROO problems have been obtained. Arvanitidis and Rosing [36] proposed a simplified equivalent reservoir representation of a multi-reservoir hydroelectric system, and Brandão [37] verified the validity of this method for modeling the multi-reservoir hydroelectric system optimal operation. Zhao evaluated the marginal utility principle of the potential energy of dammed water into hydropower depends on both reservoir storage and release [38], and was used to speed up the efficiency of solving ROO [5,39]. These studies show that the analysis of the characteristics of the problem improves the stability and high efficiency of the solving method for ROO. Water level and inflow are as key state variables (sometimes water level is also regarded as a control variable) and inputs in ROO. Therefore, the analysis of the relationship of optimal control water level and inflow is conducive to understanding and solving ROO under deterministic inflow conditions. Fuzzy cognitive map (FCM) is a graph model that uses a weighted directed graph to visualize expert knowledge [40]. In FCM, the concepts or elements are treated as nodes, and the strength of causality is expressed by the signed weights between nodes. Traditionally, domain experts allocate initial weights. Besides, Hebbian-based learning algorithms [41,42] and evolutionary algorithms [43][44][45] can improve the accuracy of the weights, and these new methods can effectively coordinate conflicts among experts. The obtained FCM model can be used to study the characteristics of complex systems. Because of these advantages, FCM has been applied to study a wide variety of complex systems, such as engineering control systems [46,47], education [48] and ecosystems [49,50].
In this paper, a new FCM with offset (FCM-O) is proposed to try to infer the causal relationship between optimal control water level and inflow of ROO. Major contributions are outlined as follows: (1) FCM-O is proposed to overcome the causal inference error caused by non-linear mapping of the activation function. In FCM-O, the activation function is not used, and the offset is introduced to better train directed weighted graphs to illustrate the specific relationship between any pair of elements. (2) The FCM-O of ROO for the Three Gorges Reservoir (TGR) is established. The causal relationships between optimal control water level and inflow are inferred using FCM-O, and they are presented as intuitive graphical forms. In addition, some relevant conclusions are obtained.
The remainder of this paper is organized as follows: Section 2 introduces the details of ROO. In Section 3, the process of obtaining the optimal control water level of ROO using dynamic programming is presented. Section 4 introduces the details of FCM and FCM-O. In Section 5, FCM and FCM-O are applied to ROO for TGR, and the experimental results are analyzed. Finally, conclusions are summarized in Section 5.

Problem Formulation
This section introduces the details of ROO, mainly including the mathematical model and several constraints.

Objective Function
The primary objective of ROO is to maximize the comprehensive benefit of water resources utilization.
The maximization of power generation is frequently selected as the optimal operation objective. The power generation represents direct economic benefits, which is a commonly used optimization objective in ROO [7,32,51]. It can be expressed as follows.
where T is the number of periods; A is output coefficient of TGR; ∆t shows interval of scheduling term; N t , H t and Q t denote output, pure water head and generating discharge in t period, respectively.

Constraints
The following constraints of ROO should be considered.
where V t is reservoir storage at t period; I t is inflow at t period; O t , Q t and S t stands for outflow, generating discharge and deserted outflow, respectively.
where H t stands for the water head; Z t is the upstream water level at t period; Z down t is tail water level described in Equation (6); The function SDR() represents the hydraulic connection between Z down t and O t .
where Z min t and Z max t are the minimum and maximum water level limits; and ∆Z is the maximum amplitude of water level variation.
where N max t (H t ) represents the maximum output, which is a function of water head; N min t is the minimum output limit.
where O min t and O max t represent respectively the minimum and maximum outflow limit.
where Z begin and Z end are initial water level and terminal water level of hydropower station, respectively.

Obtaining the Optimal Control Water Level of Reservoir Optimal Operation Using Dynamic Programming
Dynamic programming (DP) is one of the most well-known and effective methods to handle the optimization problem with the multi-stage sequential decision, which was first introduced by Bellman in 1962 [52]. The multiple-period ROO can be formulated as a multi-stage sequential decision problem, and it is extremely effective for DP to deal with nonlinear objective functions and reservoir operation constraints [53][54][55]. According to the principle of optimality, the recursive function for reservoir operation in DP can be expressed as below.
where O t,j represents the optimal control decision. Considering the water balance O t = V t−1 −V t ∆t + I t in Equation (3), the decision can be expressed by the state V t−1 and V t , the Equation (12) can be expressed as follows. Backtracking where R t (V t−1,i ) represents the benefit of the remaining period from t to T for state V t−1,i , and the reservoir capacity is V t−1,i at the beginning of period t; B t (V t−1,i , V t,j , I t ) stands for single-period benefit function. In Equation (12), the benefit of R t (V t−1,i ) and B t (V t−1,i , V t,j , I t ) are specified as many scheduling objectives, such as power generation, flood control, water supply, . . . , etc. Backtracking(V t−1,i , t) = V t,j represents the backtracking relationship between V t−1,i and V t,j . The pseudocode of DP is shown in Algorithm 1.

Input:
1: set V begin and V end ; select inflow series {I t }.

Initialization:
1: the states (reservoir capacity) are discretized 2: generate discrete set of states {V t, j }

Calculation:
1: for t = T to 1 2: select optimal decision V t, j from {V t,j } to obtain the optimal R t (V t−1,i ) 4: save the backtracking relationship Backtracking(V t−1,i , t) = V t, j 5: end for 6: end for Output: the optimal benefit R 1 (V begin ) and optimal state (decision) process {V t } Algorithm 1 and Figure 1 show the decision-making process of DP for reservoir operation. V begin and V end of TGR are set as fixed values, which are usually specified by the dispatcher. At the first and last period, there is only one state different from other periods that the reservoir capacity V is discretized into n values: {V t, j }. At t period, V t,j will be selected from {V t,j } to obtain the optimal R t (V t−1,i ) for each V t−1,i in {V t−1,i }, and save the backtracking relationship between V t−1,i and V t,j After recursive computation from T − 1 to 1 period, the optimal benefit R 1 (V begin ) and optimal state process {V * t, j } can be obtained based on the backtracking relationship previously preserved. The reservoir capacity is selected as the state in the process of DP for ROO, and water level is often used as a control attribute in practical engineering applications. The water level and reservoir capacity can be queried by the relationship between water level and reservoir capacity as Equation (15), where Czv() and Cvz() respectively represent the functional relationship between reservoir capacity and water level.

Fuzzy Cognitive Map with Offset
This section first introduces some basic concepts related to the FCM. Then a simplified but practically-improved FCM with offset will be introduced in detail.

Fuzzy Cognitive Map
FCM was proposed as a generalized cognitive map by Kosko [40]. An FCM is a graph with N nodes. Each element of the complex system under study is treated as a node. The value of node n is denoted as n C ( ) 1, 2, , n N =  . n C ranges (0, 1), and it reflects the status of the corresponding element. The directed and weighted edges connect the related and paired nodes. The direction of the edge represents the causal relationship between the nodes. An example of FCM is shown in Figure 2. The weight of the edge from node i to node j is denoted as positive ij w represents an excitatory relation from node i to node j, i.e., an increase (decrease) in i C will cause an increase (decrease) in j C , while a negative ij w represents an inhibitory relation. In addition, there is no relationship from node i to node j when the weight ij w is 0.
In summary, the optimal state (decision) process {Z * t } or {V * t,j } depends upon the inflow sequence (input), the initial and terminal water level of reservoir operation (boundary condition) and the constraints.
It shows that under the given boundary conditions and constraints, there exists a causal relationship between the optimal control water level process {Z * t } and the inflow sequence {I t }.

Fuzzy Cognitive Map with Offset
This section first introduces some basic concepts related to the FCM. Then a simplified but practically-improved FCM with offset will be introduced in detail.

Fuzzy Cognitive Map
FCM was proposed as a generalized cognitive map by Kosko [40]. An FCM is a graph with N nodes. Each element of the complex system under study is treated as a node. The value of node n is denoted as C n (n = 1, 2, . . . , N). C n ranges (0, 1), and it reflects the status of the corresponding element. The directed and weighted edges connect the related and paired nodes.
The direction of the edge represents the causal relationship between the nodes. An example of FCM is shown in Figure 2. The weight of the edge from node i to node j is denoted as w ij ∈ (−1, 1). A positive w ij represents an excitatory relation from node i to node j, i.e., an increase (decrease) in C i will cause an increase (decrease) in C j , while a negative w ij represents an inhibitory relation. In addition, there is no relationship from node i to node j when the weight w ij is 0. FCM with N variables (elements) can be expressed as a weighted matrix which is shown in Equation (16) [50,56]. According to the dataset used to construct the FCM, the value of the ith element, which is affected by the values of elements that have nonzero connections to the ith node, can be computed with Equation (17), C t are the value of the ith element at the (t + 1)th period, and the value of the jth element at the tth period, respectively. () f represents the activation function that restricts the values of the node within the range of (0, 1). There are three widely used activation functions: The signum function, the trivalent function and the sigmoid function [57]. The sigmoid function has been used by most of the studies on FCM [58]. The sigmoid activation function is defined as follows: where 0 λ is the parameter that determines the steepness of the sigmoid function at values around 0. Figure 3 shows the function graph of ( ) f x with different 0 λ . Different 0 λ make the activation function work better for the weight matrix of the actual applications. However, it is the commonly used that 0 λ is set to 5.0 for the simulated study based on data generated using FCM in the some literature [59,60].  FCM with N variables (elements) can be expressed as a weighted matrix which is shown in Equation (16) [50,56]. According to the dataset used to construct the FCM, the value of the ith element, which is affected by the values of elements that have nonzero connections to the ith node, can be computed with Equation (17), where C i (t + 1) and C j (t) are the value of the ith element at the (t + 1)th period, and the value of the jth element at the tth period, respectively. f () represents the activation function that restricts the values of the node within the range of (0, 1). There are three widely used activation functions: The signum function, the trivalent function and the sigmoid function [57]. The sigmoid function has been used by most of the studies on FCM [58]. The sigmoid activation function is defined as follows: where λ 0 is the parameter that determines the steepness of the sigmoid function at values around 0. Figure 3 shows the function graph of f (x) with different λ 0 . Different λ 0 make the activation function work better for the weight matrix of the actual applications. However, it is the commonly used that λ 0 is set to 5.0 for the simulated study based on data generated using FCM in the some literature [59,60]. FCM with N variables (elements) can be expressed as a weighted matrix which is shown in Equation (16) [50,56]. According to the dataset used to construct the FCM, the value of the ith element, which is affected by the values of elements that have nonzero connections to the ith node, can be computed with Equation (17), where ( 1) i C t + and ( ) j C t are the value of the ith element at the (t + 1)th period, and the value of the jth element at the tth period, respectively. () f represents the activation function that restricts the values of the node within the range of (0, 1). There are three widely used activation functions: The signum function, the trivalent function and the sigmoid function [57]. The sigmoid function has been used by most of the studies on FCM [58]. The sigmoid activation function is defined as follows: where 0 λ is the parameter that determines the steepness of the sigmoid function at values around 0. Figure 3 shows the function graph of ( ) f x with different 0 λ . Different 0 λ make the activation function work better for the weight matrix of the actual applications. However, it is the commonly used that 0 λ is set to 5.0 for the simulated study based on data generated using FCM in the some literature [59,60].

Fuzzy Cognitive Map with Offset
The activation function f (x) can restrict the node values within the range of (0, 1). But the effect of f (x) is greatly affected by the parameter λ 0 , and f (x) converts N i w ij C j (t) nonlinearly into (0, 1), which may cause data distortion. FCM with an offset is proposed to try to overcome the above problems. In FCM-O, the activation function is not used, and the offset is introduced to better learn the structure of a directed weighted graph to illustrate the concrete relationships between any pair of elements. The weighted matrix of FCM-O is the same as that of FCM as Equation (16). The dynamics of the node values in the FCM-O is simulated by using the following equation: where O f f set i is the offset of the ith element. The application of offset is to overcome the causal inference error caused by non-linear mapping of activation function could improve the fitting accuracy of the FCM, and the result of causal inference will be more intuitive without using the activation function. The processes of learning the structure of FCM and FCM-O is the same, which will be described in detail in the following Section 4.3.

Algorithm for Learning the Structure of FCM: Differential Evolution Algorithm
Most of the studies apply data-driven FCM learning algorithms focused on minimizing the simulation error [61][62][63] between an output sequence and historical data. A large number of meta-heuristic FCM learning algorithms exist [63], and can be found in the literature about differential evolution (DE) algorithm [64].
DE is a simple yet practical method for optimization problems proposed by Price and Storn [65]. Due to its simplicity, easy implementation, fast convergence and robustness, the DE algorithm has attracted more and more attention, and has been widely applied successfully, such as flow shop scheduling, piezoelectric performance optimization, economic optimization design, etc. [7,66,67]. A virtual population is used to form an initial solution in feasible space, and then recursive implementations of differential mutation, crossover, and greedy selection are executed until the termination condition is met. The solutions in the population evolve into a solution suitable for the dataset in this evolutionary process.
Finally, an optimal solution that can reasonably and accurately represent the FCM is formed. The solutions in the population generally represent the weight matrix and some parameters in FCM. In the greedy selection procedure, the min data error (fitness function) is used to evaluate the solutions in the population. The objective function in this research is shown as Equation (20).
where T is the period number of datasets. C n (t) and C n (t) are the value of the nth element at the tth moment in the dataset and the corresponding value returned from FCM, respectively. The smaller the data error is, the better.
LSHADE [68] contains the external archive and the mutation strategy 'DE/current-to-pbest/1' like JADE [69] and historical memory of successful control parameter settings like SHADE [70], and the strategy of linear population size reduction is applied. It secures the first rank (see http: //www.ntu.edu.sg/home/EPNSugan/index_files/CEC2014/-CEC2014.htm) on real-parameter single objective numerical optimization in CEC2014 competition [71]. Due to the excellent performance of LSHADE in numeric optimization problems, it is selected to learn the structure of FCM. The parameter setting and specific implementation of LSHADE are referred to Ref [68].
Mutation is a change or perturbation with a random element on the population to locate better solutions. DE creates a mutant vector → v i,G by utilizing one or more difference vectors. The mutation strategy in original DE is "rand/1", which is expressed in Equation (21).
Other common DE mutation strategies are as follows: "rand/2": "best/1": "best/2": "current to best/1": where the indices r1 − r5 represent the random and mutually different integers generated within the range {1, NP} and different from the index i.

→
x best,G is the best individual in a current generation. Each strategy has a different ability to maintain the diversity of the population, which may increase/reduce the rate of convergence in the process of evolution.
LSHADE employs the mutation strategy 'DE/current-to-pbest/1' proposed in JADE, which is shown in Equation (26). 'DE/current-to-pbest/1' utilizes an optional external archive to diversify mutant vectors, and this achieves good performance over many benchmark functions. → x pbest,G in Equation (26) denotes certain randomly-chosen vectors from the top 100 × p% individuals. DE/current-to-pbest/1 depends on the control parameter p in order to balance exploitation and exploration (small p behaves more greedily), and p is a constant set to 0.05 in JADE. JADE defines a set of archived inferior solutions A, and current population solutions P. → x r2,G in mutation strategy is selected from the union P∪A. The archive A is initialized empty, and then inferior solutions (individuals) are added to the archive in each generation. The upper bound of archive size is fixed, and if the size exceeds the fixed threshold, these redundant solutions are randomly selected and removed from the external archive A to keep the upper bound of archive size equaling to the fixed threshold.
After mutation, LSHADE generates a trial vector → u i,G by a crossover operation to enhance the potential diversity of the population. In the basic version, the binomial crossover is a commonly used crossover operator operation in DE. It crosses the parent individual where j rand is an index of uniformly randomly chosen from [1, D], which ensures that the trial vector The corresponding solution and control parameters are called successful solution and successful control parameters, respectively. Otherwise, they are called unsuccessful update.
For the linear population size reduction, formula (29) shows the dynamic change of population size in LSHADE. NP min is set to the smallest possible value, such that the evolutionary operators can be applied in the case of L-SHADE, NP min = 4 because the DE/current-to-pbest/1 requires 4 individuals in Equation (26). NFE is the current number of fitness evaluations, and MAXNFE is the maximum number of fitness evaluations. Whenever NP G+1 <NP G , the (NP G − NP G+1 ) worst-ranking individuals are deleted from the population, and the external archive size |A| changes adaptively according to the dynamic population size, |A G | = round(NP × r arc ), r arc is the archive factor. Whenever the size of the archive exceeds the predefined archive size |A|, the randomly-selected elements are deleted to make space for the newly-inserted elements.
Control parameters in LSHADE are also updated in an adaptive manner as SHADE. F obeys a Cauchy distribution in which the location parameter is denoted by µF and the scale parameter equals to 0.1, F ∼ C (MF, 0.1). Cr obeys a normal distribution in which the mean value is denoted by MCr and standard deviation equals to 0.1, Cr ∼ N (MCr, 0.1). If the values generated for F and Cr are outside the range (0, 1), Cr will be truncated to (0, 1), and F will be truncated to be 1 if F ≥ 1, or regenerated if F ≤ 0. The updating schemes of control parameters F and Cr are listed in Equations (30) and (31), respectively.
Success values of F and Cr in each generation are recorded and used to update SF and SCr. The updating schemes of control parameters MF and MCr are listed in Equations (32) and (33), respectively. As MCr is updated, if MCr h,G = "⊥" (where "⊥" denotes a special, "terminal value") or max(SCr) = 0 (i.e., all elements of SCr are 0), MCr h,G+1 is set to "⊥". Thus, if MCr is assigned the terminal value "⊥", then MCr will remain fixed at "⊥" until the end of the search.
This has the effect of locking Cr to 0 until the end of the search, causing the algorithm to enforce a "change-one-parameter-at-a-time" policy [72], which tends to slow down convergence, and is effective on multimodal problems.
where index h (1 < h < H) determines the position in the memory to update. The index h is initialized to 1 at the beginning. h is incremented whenever the entry (MF h , MCR h ) is updated once. If h > H, h is set to 1. In the update Equations (32) and (33), the weighted mean mean WA (SCr) is computed according to Equation (35) referred to by Peng et al. [73]. The weighted Lehmer mean mean WL (SF) is computed using the formula below, and as with mean WA (SCr):

Description of Research Area
The Yangtze River is the third longest river in the world, the largest river in China. The Three Gorges Reservoir (TGR) is the key backbone of the Yangtze River Basin (See as Figure 4), and the operation and management of TGR will inevitably affect the river basin ecosystem of the upper and lower reaches. The main parameters of TGR are listed in Table 1. In order to defend the coming flood in the flood season, the TGR empties storage capacity, and its water level dropped to 145 m at June 10. Then TGR begins to store water in September 1.

Description of Research Area
The Yangtze River is the third longest river in the world, the largest river in China. The Three Gorges Reservoir (TGR) is the key backbone of the Yangtze River Basin (See as Figure 4), and the operation and management of TGR will inevitably affect the river basin ecosystem of the upper and lower reaches. The main parameters of TGR are listed in Table 1. In order to defend the coming flood in the flood season, the TGR empties storage capacity, and its water level dropped to 145 m at June 10. Then TGR begins to store water in September 1.

Dataset Acquisition and Preprocessing
The main purpose of the study is to infer the causal relationship between the optimal control water level process { * t Z } and the inflow sequence { t I } in ROO. Section 3 describes in detail how to obtain an optimal control water level of reservoir operation by DP. The time step was monthly.

Dataset Acquisition and Preprocessing
The main purpose of the study is to infer the causal relationship between the optimal control water level process {Z * t } and the inflow sequence {I t } in ROO. Section 3 describes in detail how to obtain an optimal control water level of reservoir operation by DP. The time step was monthly. Therefore, the scheduling period begins at June 1 and ends at June 1 of the next year. The initial water level Z begin i and terminal water level Z end i are set to 145 m. There is the historical inflow data of TGR for 55 years from 1959 to 2014. The objective of reservoir optimal operation is to maximize power generation as Equation (1). Figure 5 shows the inflow sequence {I t } and the optimal control water level process {Z * t } obtained by DP. From June to September, the water level of TGR is limited to the flood limit level (145 m) to prevent floods. Therefore, our research period focuses on the non-flood season, that is, September to May of the next year. {Z * t } and {I t } are illustrated in Figure 6 under the condition of a certain inflow sequence. Z f represents the flood limit level. Therefore, the scheduling period begins at June 1 and ends at June 1 of the next year. The initial water level begin i Z and terminal water level end i Z are set to 145 m. There is the historical inflow data of TGR for 55 years from 1959 to 2014. The objective of reservoir optimal operation is to maximize power generation as Equation (1). Figure 5 shows the inflow sequence { t I } and the optimal control water level process { * t Z } obtained by DP. From June to September, the water level of TGR is limited to the flood limit level (145 m) to prevent floods. Therefore, our research period focuses on the non-flood season, that is, September to May of the next year. { * t Z } and { t I } are illustrated in Figure 6 under the condition of a certain inflow sequence. f Z represents the flood limit level.
According to the regular "Optimal Operation Scheme of the Three Gorges" drawn up by Chengdu Engineering Corporation Limited, the water level control generally follows these regulations: (1) For power generation, the water level of TGR should be higher than 145 m, which is the dead water level. In addition, TGR should keep lower than the normal water level 175 m. (2) From July to early September, the TGR runs according to the flood control mode.
(3) TGR begins to store water in September, and reaches 175 m by late October. TGR had better fill up quickly to improve the efficiency of power generation.       According to the regular "Optimal Operation Scheme of the Three Gorges" drawn up by Chengdu Engineering Corporation Limited, the water level control generally follows these regulations: (1) For power generation, the water level of TGR should be higher than 145 m, which is the dead water level. In addition, TGR should keep lower than the normal water level 175 m. (2) From July to early September, the TGR runs according to the flood control mode.
(3) TGR begins to store water in September, and reaches 175 m by late October. TGR had better fill up quickly to improve the efficiency of power generation.
After the inflow sequence {I t } and the optimal control water level process {Z * t } for 55 years from 1959 to 2014 obtained by DP, all the data of {Z * t } and {I t } were preprocessed with Equations (37) and (38) to normalize the influence of each element and avoid the dominant elements whose absolute values are much larger than those of the other attributes. Z t,y and Z * t,y (t = t 1 , . . . , t T−1 ; y = 1, . . . Y) are the normalized and pre-normalized optimal control water level of tth period with the yth inflow sequence. Z min and Z max are the minimum and maximum values of the optimal control water level under all the different inflow sequence conditions, respectively. T is the number of research periods. Y is the number of inflow sequence. I t,y and I t,y (t = t 1 , . . . , t T ; y = 1, . . . Y) are the normalized and pre-normalized inflow of tth period with the yth inflow sequence. I min and I max are the minimum and maximum values of the optimal control water level under all the different inflow sequence conditions, respectively. The month corresponding to t = t 1 , . . . , t T−1 is shown in Table 2. All the data were randomly divided into two parts (3/4 and 1/4) for training and testing. The final result is the FCM or FCM-O model with minimum data error.

Case Study and Discussion
In this case, FCM and FCM-O are established to infer the causal relationship between optimal control water level {Z * t,y } and inflow {I t,y } of ROO for TGR. In order to simplify the model, the relationship between inflows in each month is neglected. In addition, the recursive formula of DP in Section 3 shows that {Z * t,y } mainly depends on {I t,y } under given boundary conditions and constraints. Therefore, the model description of the above FCM and FCM-O can be simplified to the following expressions. (40) where Z tZ,y is the tZth optimal control water level at the yth inflow sequence returned from FCM or FCM-O. w tZtI represents the weight of the edge from the tIth inflow to the tZth optimal control water level. The corresponding objective function in this case is shown as Equation (25).
LSHADE is selected to learn the structure of FCM, and the number of max fitness evaluations is set to 10,000 × D. D is the dimension of the problem. From Figure 6, it can be seen that the number of node Z * t,y is 8, and the number of node I t,y is 9 (T = 9). So D = 8 × 9 + 8 = 81 in FCM and FCM-O (λ 0 of f (x) in FCM and O f f set tZ in FCM-O are needed to search for each Z * t,y ). Table 3 shows the min data error of FCM and FCM-O for training and testing. Compared with FCM in the min data error, FCM-O reduces 11.11% and 7.14% in the training and the testing, respectively. It shows that FCM-O is superior to FCM in terms of its minimum data error. The weight matrix and parameter of training results of FCM and FCM-O are shown in Tables 4 and 5.  It can be seen from Figure 5a that Z * 9 , Z * 10 and Z * 11 maintain the same water level under different inflow sequences. TGR begins to store water in early September. According to the regulations, the water level at the end of September is not higher than 162 m, while Z * 9 reaches the upper boundaries of the constraint at 162 m. TGR is generally full to 175 m of normal water level from October to December, while Z * 10 and Z * 11 reaches the upper boundaries of the constraint at 175 m. It illustrates that TGR should impound water as fast as possible during the impoundment period, raising the head of power generation so as to increase the total power generation, which is consistent with the conclusions of many literatures [5,32,38,74]. Figure 5b shows that I 9 ranges in (12,000, 40,000) m 3 /s, while I 10 ranges in (11,000, 28,000) m 3 /s and I 11 ranges in (8,000, 16,000) m 3 /s. These indicate the value of Z * 9 , Z * 10 and Z * 11 are more affected by upper boundaries of the water level constraint constraints than by inflow. However, the weighting coefficients of I t,y to Z * 9 , Z * 10 , and Z * 11 in FCM are 1 or −1 in Table 4, which shows that Z * 9 , Z * 10 and Z * 11 are greatly affected by I t,y . The conclusions obtained by FCM do not match the accepted conclusions that Z * 9 , Z * 10 and Z * 11 are more affected by constraints than by inflow, and the weighting coefficients of I t,y to Z * 9 , Z * 10 and Z * 11 in FCM-O are 0, which seems reasonable. The reason for the above causal inference errors in FCM lies in the inability of the activation function f (x) to deal with variables with offsets that are not 1 or −1. All of those shows that FCM-O is more competent than FCM in the causal relationship between optimal control water level and inflow in ROO both in min data error and reasonableness of results.
To analyze the results of FCM-O, the weight matrices corresponding to all subfigures are shown in Figure 6. Because the relationship between inflows is neglected, and the optimal control water level depends on the sequence of inflow, and are unrelated to the adjacent water level, the optimal control water level only pays attention to the in degree and offset, and all out degrees are 0. Table 6 shows the in degree and offset of each Z * t of Figure 7. Z * 9 , Z * 10 , Z * 11 and Z * 12 need to reach higher operating water levels, which are more affected by constraints than by inflow. So there is no inflow but only offset acting on Z * 9 , Z * 10 , Z * 11 and Z * 12 , as shown in Figure 7a-d. Z * 1 , Z * 2 and Z * 3 reflect the optimal control water level in the early stage of the water supply period, and Z * 1 , Z * 2 and Z * 3 range in (145, 175) m. In the early stage of the water supply period, TRG should reduce the discharge as much as possible and maintain a higher head, so as to increase power generation under the premise of meeting the minimum discharge constraint. Therefore, Z * 1 , Z * 2 and Z * 3 are mainly affected by the inflow in recent months. The larger the inflow, the smaller the amount of water that the reservoir needs to release, and correspondingly, Z * 1 , Z * 2 and Z * 3 will be higher. The inflow in recent months has shown positive effects to Z * 1 , Z * 2 and Z * 3 . Z * 4 reflects the optimal control water level at the end of the water supply period (before the flood season). There are many elements to be considered in Z * 4 . First, Z * 4 needs to maintain a high water level to ensure a higher head, thereby increasing power generation. Secondly, Z * 4 needs to avoid too high a water level to produce discarded water when emptying the reservoir before flood season. Therefore, Z * 4 has a narrow range and is less affected by inflow in different months as shown in Figures 5a and 6h. Z . In addition, we can use this relationship to estimate the optimal control water level with inflow The results in this research are consistent with some abovementioned literatures. Moreover, our studies made these conclusions quantitative. Meanwhile, we found some details that were overlooked previously, such as the conclusion that Z * 9 , Z * 10 and Z * 11 are more affected by constraints than by inflow, and the inflow in recent months has shown positive effect to Z * 1 , Z * 2 and Z * 3 . In addition, we can use this relationship to estimate the optimal control water level with inflow in FCM-O.

Conclusions
A new FCM with offset was proposed to try to infer the causal relationship between the optimal control water level and the inflow of ROO. To test its performance, the FCM and FCM-O of ROO for TGR are established. The optimal control water levels of TGR were obtained by DP with the historical inflow data for 55 years from 1959 to 2014. Then all the data of this optimal control water level and inflow were randomly divided into two parts (3/4 and 1/4) for training and testing. The FCM or FCM-O model were trained with minimum data error by LSHADE. The experimental results show that, compared with FCM in the min data error, FCM-O reduces 11.11% and 7.14% in the training and the testing, respectively. In addition, FCM-O can deduce that the optimal control water level from September to December needs to be raised as much as possible, which are mainly restricted by constraints and are not affected by inflow. The optimal control water level from January to March also needs to maintain a high water level, and the inflow in recent months has a positive effect upon it. The optimal controlled water level at the end of April is at the end of the water supply period and near the flood season. The optimal control water level at the end of April needs to maintain an appropriate high water level, so as to avoid the risk of discarding water while guaranteeing high water head to increase power generation. In addition, the relationship obtained by FCM-O can be used to estimate the optimal control water level with inflow.