Simulation and optimization of control system operation and surface water allocation based on system dynamics modeling

Iran is in the dry belt of the world and its rainfall is one-third of the average global rainfall. Drought periods in Iran are frequent based on the records. Therefore, proper planning such as taking effective steps for the optimal use of water is required to minimize the damage of water shortage. The system dynamics (SD) method is one of the most powerful management tools which has been used in recent years to simulate complex water resource systems. Hence, this method is applied in the present study to simulate the performance of Dez Reservoir. It is important to note that Dez Reservoir plays a significant role in supplying water for 125,000 ha downstream irrigation and drainage network along with its effective role in providing hydroelectric power and flood control.

The system dynamics in the decision-making process of the reservoir management not only have the ability to analyze the complex systems based on reality, but it also provides the opportunity for the user's participation in the development of the model. Ahmad & Simonovic (2000) developed a model for flood management in the Shellmouth Reservoir in Canada using the SD method. Hassanzadeh et al. (2012) evaluated the main factors in Urmia Lake's level reduction using the SD model. After validating this model, the results showed that changes in inflows due to climate change and overuse of surface resources were the main factors in reducing the surface area of the Urmia Lake by up to 65%. The construction of four different dams and reduction of precipitation resulted in decreasing the lake level up to 25 and 10%, respectively. Madani & Mariño (2009) used the SD method to manage the Zayandehrood River basin. Results of the model for different scenarios showed that inter-basin transfer is not the only solution to the river deficits. Various alternatives such as demand management and population control can be used along with inter-basin water transfer to solve the basin water crisis more effectively by increasing the water storage capacity and controlling the groundwater discharge.

Kotir et al. (2016) studied water resource management in the Volta River Basin in Ghana by designing three scenarios of infrastructure development, wheat crop expansion, and arid conditions using the SD method. This study indicated that the scenario of infrastructure development is an important step towards improving the region's agriculture as well as conserving water resources. Xi & Poh (2013) used the SD models to analyze Singapore's integrated water resources and demonstrated that the system dynamics is a suitable decision support tool for sustainable water resource management. Moreover, Nozari et al. (2014) simulated the right part of the drainage network of Abshar in Isfahan as well as the cropping pattern in the concerned area using the SD method. The results of the model validation showed that the model has acceptable accuracy in simulating an irrigation network, cropping pattern, and other defined scenarios. Several recent studies have used the SD method to develop dynamic systems and simulation models in different river basins or watersheds around the world (Niazi et al. 2014; Kotir et al. 2017). Li et al. (2018) simulated water resource planning in Zhengzhou, China, using the SD model. The model consisted of four subsystems of population, economic, water demand, and water supply. The difference between supply and demand in Zhengzhou is very conspicuous. The results reveal that in 2020 the demand for water will reach 2.57 billion cubic meters in Zhengzhou.

Furthermore, reservoir operation optimization is studied using various algorithm techniques by many researchers. For example, heuristic approaches and evolutionary algorithms have been studied by Zhang et al. (2016), Niu et al. (2018), Fang & Popole (2019), Moeini & Afshar (2019) and Azad et al. (2020). Moeini & Afshar (2019) presented two hybrid ant colony optimization algorithm methods (ACOA-NLP and ACOA-TGA-NLP) for the layout and pipe size optimization of sewer network. In this study, the computational effort is significantly reduced in comparison with alternative methods. Dehghani et al. (2019) applied the Grey wolf optimization (GWO) method coupled with an adaptive neuro-fuzzy inference system (ANFIS) to forecast the hydropower generation. Azamathulla et al. (2008) applied an integrated linear programming (LP) and genetic algorithm (GA) model to optimize the allocation of the available water in the Chiller reservoir system in India. They concluded that GA algorithm performs much better compared to the LP model. Fang & Popole (2019) utilized the developed PSO algorithm for multi-objective optimization of the Shui-Kou hydroelectric reservoir in Minjiang, China. This research showed that the model can achieve the optimum schedule by considering environmental and electricity generation benefits; additionally, an optimal operation plan can improve environmental protection by 6.56% and reduce electricity production costs by 0.47%. Niu et al. (2018) used the multi-objective PSO algorithm in order to optimize the exploitation of the Lancang hydroelectric reservoir in the southwest of China. The results indicated that the multi-objective PSO algorithm can provide satisfactory planning results for the reservoirs and various other cases. Zhang et al. (2016) investigated the optimal exploitation of the multi-reservoir hydraulic systems using PSO optimization and their results showed improvement of the performance and operation of this algorithm in high-quality volume storage.

In the modeling process, main sources of uncertainty are in input data, model structure and parameter values. Hence, it is important to quantify the uncertainty in the model outputs. Fan et al. (2019) developed a full-subsampling factorial copula and a full-subsampling factorial analysis approach to quantify parameter uncertainties and to diminish the effect of the sample size, respectively. Li et al. (2019) proposed a fuzzy gradient chance-constrained evacuation model to address different uncertainties in evacuation management and planning. One of the most widely used methods for calibration and uncertainty analysis of hydrological model parameters is Generalized Likelihood Uncertainty Estimation (GLUE). The reason for using this method is its simplicity, ease of performance, and less strict statistical assumptions about model errors (Cho et al. 2019).

Studying the conducted research shows that the dam reservoir modeling and its performance simulation using the SD method has acceptable accuracy. However, so far, the performance of the reservoir, the regulatory dam, and the diversion dam of the surface water control system has not been simulated simultaneously using this method. Therefore, in the present study, the performance of the SD method in the management and simulation of the productivity of the water resource system for the storage reservoir, regulatory dam, and diversion dam will be investigated simultaneously. Also, the accuracy of this method in the management of the mentioned integrated system will be compared with that of the meta-heuristic particle swarm optimization (PSO) algorithm.

The Dez basin is part of the second-grade basin of the great Karun. There are 27 water resource development projects in this basin, one of which is the Dez Reservoir dam. The minimum and maximum storage capacity of the dam are 466.73 and 3,460 million cubic meters, respectively (Adib et al. 2012). The main purpose of the construction of this dam was to control floods, and to regulate water for irrigation and agricultural purposes. In order to control the fluctuations of the outflow water from the Dez Reservoir dam, the Dezful regulatory dam was constructed at 32°24′33″ N and 48°24′56″ E and it was located 5.4 km to the north of the Dezful old bridge to control the outflow water from the Dez Reservoir dam and to supply inflow water to the Dez irrigation and drainage network based on 200 cubic meters per second. Sabili canal, with a capacity of 13 cubic meters per second, was constructed 2 km upstream of the dam. Also, in order to divert the Dez River water to the main west (with the capacity of 157 cubic meters per second) and east canals (with the capacity of 93 cubic meters per second) gravitationally, supply the downstream water, and cross 6,000 cubic meters per second of flood flow, the diversion dam was constructed. The dam, with a length of 394 m and a height of 4 m is located 6 km away from the regulatory dam's downstream (Figure 1). The characteristics of the reservoir dam, regulatory dam, and diversion dam are presented in Table 1.

In the above equation, a and b are the constant coefficients obtained from the area-storage curve of the reservoir which are 8.7054 and 0.0207, respectively. The monthly evaporation heights of the Dez dam station and the Dezful regulatory dam are presented by evap term. S_t and S_t+1, and are respectively the reservoir volume at the beginning and the end of the period.

The required flow to preserve the habitat of various species of aquatic organisms and plants to conserve their life cycle is called the minimum environmental demand. After exploitation, each dam should meet its downstream environmental demands. The Montana method was used for the quantitative estimation of the Dez dam's environmental demands. In this method, by considering the aquatic life conditions, the minimum amount of water release is calculated as a specific percentage of the annual average river discharge (Table 3). In order to estimate the environmental demand, in the Montana method, 10% of the mean monthly discharge for the rainy months and 30% of the mean monthly discharge for the low-water months of the year were used as the long-term statistical series.

The SD method is a method which is based on systematic thinking. It is used for studying and improving learning in complex systems. This approach is fundamentally feedback-based. In this method, four tools of stock, flow, interfaces, and converters are used (Niazi et al. 2014). In a rapidly changing water resource system, dynamic models can be useful analytical tools for studying the behavior of the variables through time. The SD, in addition to being able to describe complex systems based on reality, enables the user's participation in developing the model so that the best decision can be made.

In the present study, the Ventana system (Vensim) is used to evaluate system dynamics. This tool is an object-oriented modeling environment that makes it possible to create complex models more easily compared to conventional programming languages. Using this software, the Causal Loop Diagram (CLDS) can be converted to Stock and Flow Diagrams (SFDS) (Kotir et al. 2016).

SD modeling is not able to provide a frame or specific methodology for the design and diagnosis of organizational structures in the sense of mutual relation among organizational actors. Absorption of variety by an organization is another limitation of SD modeling. It is important to note that variety is the number of states a system can assume and it has been used as technical term for complexity (Schwaninger 2020). Another significant limitation for SD models is lack of enough skill and experience to build the desired model. In other words, learning how to build models and examine them not only needs years of training, but also a system dynamic collaborator is required to help the community psychologists to apply necessary techniques.

Based on the SD modelling, we are able to model a problem in different scales by focusing on the identification of main stock players at a specified resolution. However, we cannot ignore the influence of parameters and auxiliary variables. Although thinking and modeling at different scales is an important part of this approach, there will not be a formal procedure to cope with the external complexity. Furthermore, SD modelling generates more realistic representations of complex systems in comparison to the path diagrammatic models and traditional social modelling approaches, but they cannot capture the full complexity of a process. Hence, it is important to pick appropriate elements and interactions that help modelers to create a more useful insight to define a complex problem into a single model. Although system dynamics has the potential to simulate a complex system, it cannot do so perfectly due to lack of data availability, systems understanding, and systems uncertainties.

Although a main reason for using these types of models is to make the assumptions as reliable as possible, the SD and other types of models still cannot completely satisfy all the required critics.

In any SD modeling project, the first step is to determine the system structure, which includes positive and negative relationships among variables, feedback loops, system structure, and delays (Winz et al. 2009). Using causal loops, the affects and influences of each of the variables on one another can be indicated. The causal relationship is shown by two positive and negative signs, if the changes of the independent and the dependent variable are in one direction, the causal relationship will be positive (R), and if their changes are in opposite directions, the causal relationship will be negative (B) (Kotir et al. 2017). The causal loops of the storage reservoir, regulatory dam, and diversion dam are presented in Figure 2. For example, in the B1 loop, raising water elevation increases the reservoir surface and so the volume of evaporation will be increased. Increase in evaporation reduces the storage capacity of the reservoir and eventually decreases the height of the water behind the dam and creates a negative ring.

After drawing the causal loop graphs, it is necessary to develop flow charts. The integrated state flow diagram of the Dez Reservoir dam, regulatory dam, and diversion dam is presented in Figure 3. In this stock and flow diagram, the capacity of the reservoir, regulatory dam, and diversion dam are considered as stock and they are plotted against other variables such as evaporation, overflow, release rate (to meet the demands of the turbine), release values (for agricultural use through each of the eastern, western and Sabili canals), and the consumption related to environmental demands, as flow. The characteristics and detailed information of all desired variables used in the stock-flow diagram are presented in Table 4.

The PSO algorithm is a method for iterative-based optimization which the system considers as the initial value to generate a series of random solutions. In this solution space, particles are forced to search for a better, and ultimately optimal, solution (Hu & Li 2012). The invention of the PSO optimization algorithm was inspired by the social behavior of bees, fish, or a flock of birds; as a matter of fact, the algorithm mimics the behavior of these communities. In the PSO algorithm, each solution is one bird of the flock of birds which is called particle; in this framework, birds, along with enjoying artificial intelligence, develop some social behaviors, and the coordinate of their movement is towards a destination (Montalvo et al. 2008). The basis of the PSO algorithm is that at any given moment each particle adjusts its place in the search space according to the best place it can be located. Each particle has two characteristics: position and velocity. After the dispersion of the particles, each particle communicates the position and the value of its objective function to other members. Then, the position of the best particle is identified based on which particles are afterwards redirected to the new position at a certain speed. The main advantage of the PSO is that it is easy to implement, and it requires fewer parameters. The flowchart of the PSO algorithm is presented in Figure 4.

In these equations, R_t is the monthly release rate, D_t is the monthly demand, D_max is maximum monthly demand, S_t is the reservoir storage volume at the beginning of the t period, S_t+1 is the reservoir storage volume at the end of the t period, I_t is monthly input rate of the reservoir, E_t is monthly evaporation volume. S_max and S_min are respectively the maximum and minimum storage capacity of the reservoir during the simulation period.

In the above equations is the observed discharge, is simulation discharge, is mean observed discharge, is mean simulation discharge, and n is number of data.

As stated before, one of the main objectives in this research is the assessment of uncertainty associated with the model outputs.

The basic steps of the hired methodology are as follows:

1. Seasonal autoregressive integrated moving average (SARIMA) was used to predict the monthly inflow into the Dez dam and the monthly evaporation heights of the Dez dam station. Based on the observations, SARIMA(2,1,1)(1,0,3)₁₂ was the best model for predicting monthly inflow, and the appropriate form of the monthly evaporation time series was SARIMA(1,1,3)(1,0,4)₁₂.

2. Artificial sets of data points were generated many times by using selected model and Monte Carlo simulation.

3. The SD and PSO models was run for all different input data sets.

4. The likelihood function was calculated using Equation (8) for each dataset (Beven & Binley 1992):

   

   (8)

   where is the amount of demand at time t, is the output of the model for dataset and n is the number of demand data.
5. Model output uncertainty was determined for both models by likelihood weight (Equation (9)).

   

   (9)

   is the probability or likelihood weight of the ith dataset and N is the number of datasets. The sum of these weights is equal to one and they represent the probability density function.

6. Model output uncertainty is quantified within the 95% confidence interval. In other words, data within 2.5% upper limit and 2.5% lower limit are considered as outlier data and removed from the database.

7. The ratio of the average distance between the above levels and the standard deviation of the corresponding measured variable, is calculated by the d-factor (Equation (10)) (Abbaspour et al. 2007):

   

   (10)

   

   (11)

   where is the standard deviation of the demands. is the average distance between the upper (⁠⁠) and the lower (⁠⁠) parts of the cumulative distribution. A value of less than 1 is a desirable measure for the d-factor.

In order to validate the model results, the total downstream demand was considered equal to the actual release in the years 2005–2006 to 2013–2014 and the desired model was implemented. Figure 5 shows the release values from: (a) Dez dam reservoir, (b) regulatory dam, and (c) diversion dam against the values simulated by the SD model. For better comparison of the results, the bisector line was also drawn. It is observed that in all three cases data are close to the bisector line, and this indicates the high accuracy of the model in the simultaneous simulation of these three reservoirs.

Table 5 shows the results of calculating the statistical indices of R², RMSE, and SE for variables of the water height behind the reservoir dam, the reservoir storage volume, release and overflow from the Dez Reservoir dam. It is observed that the dimensionless SE index for these variables is 0.02, 0.18, 0.08, and 0.06, respectively, indicating the high accuracy of the model in simulating the water behind the Dez Reservoir dam.

The statistical indices of R², RMSE and SE for the outflow of the regulatory dam after supplying its upstream Sabili canal demands are 0.98, 51.68 and 0.13, respectively. In addition those indexes are 0.92, 71.38, and 0.37 for the outflow of the diversion dam after meeting the main western and eastern canal requirements and the downstream environmental demands. In the Dez diversion dam, due to its lower height, calculation of the height-volume was associated with greater error which resulted in the lower accuracy of the SD model in simulating the diversion dam compared to Dez and the regulatory dam. However, in general, the results of calculating the presented statistical indices of R², RMSE, and SE show the high accuracy of the SD model in simulating the performance of the regulatory dam and the diversion dam.

The real downstream demand was defined as an input for the model. The amount of water allocation and release from each of the reservoir, regulatory, and diversion dams was calculated using the SD model. Moreover, in order to compare the results of the model, the release rate from each of the above-mentioned dams was optimized using the PSO optimization algorithm, and its results were compared with those of the SD model. The results of the implementation of these models for each of the three dams are given below.

Table 6 shows characteristics of the fitted patterns to optimize the downstream water allocation using the PSO algorithm. According to the target function values (min TSD) and the statistical indices of R², RMSE, and SE, the chosen model is model number (5) which indicated the least discrepancy between demand and release values.

Figure 6 shows the actual demand values by releasing in the three modes as: (a) observational, (b) PSO optimization, and (c) SD simulation, for the Dez. Also, the shortage values in supplying the demands of the desired reservoir for each of the three modes in all the mentioned statistical years are presented in Table 7. It is observed that the release values by the SD model are more in line with the downstream demands, and the deficiency values in this method are much lower than those of the observational and PSO optimization mode.

Table 8 indicates the desired statistical indices’ values for observational release, PSO optimization, and SD simulation. The values of these indices indicate that the SD model more accurately simulates the distribution of existing water resources to meet downstream requirements.

As mentioned earlier, the Sabili canal is located behind the regulatory dam and its input is provided by the regulatory dam. In addition, the regulatory dam should release part of the water to meet the demands of the east and west canals, towards the Dez diversion. Figure 7 shows the values of the Sabili canal demands and the downstream release values of the regulatory dam in the three simulation models, (a) observational, (b) PSO optimization, and (c) SD, in order to be compared with the downstream demands. According to this figure, it can be said that during the statistical period of 2005–2014, there was no real release in line with the demands of the Sabili canal and the downstream of the regulatory dam, and there was a significant difference between these two parameters. However, results of the PSO optimization and SD simulation demonstrate better management in meeting the demands of the subsystem canal and the downstream regulatory dam.

Table 9 shows the deficiency values in the Sabili canal and the downstream regulatory dam for the real conditions, PSO optimization, and SD simulations for the statistical study period. Moreover, the results of the statistical indices that meet this demand are presented in Table 10. Based on these tables, it can be concluded that the SD model can play an effective role in supplying the demands of both the Sabili canal and the downstream regulatory dam.

The water released from the regulatory dam is accumulated in the diversion dam and it is directed to the west and east canals. Figure 8 shows the amount of water directed to the west and east canals for the three modes of: (a) observational, (b) PSO optimization, and (c) SD simulation. According to this figure, it is observed that during the statistical period of 2005–2014, the actual release did not supply the demands of the two mentioned canals and there is a significant difference between them. Comparing the results of the PSO optimization and SD simulation, it can be stated that the results of the SD model are more in line with the downstream demands.

The deficiency values in supplying the demands of the west and east canals for the three modes of the real conditions, PSO optimization, and SD method through the considered statistical period are presented in Table 11. The corresponding statistical indices for meeting the demand are presented in Table 12. According to the results, it can be said that the SD method is more capable in surface water management.

The results showed that the SD-simulated releases met downstream demands and the SD method has high accuracy in simulating the performance of the surface tank system (Figure 9 and Table 13).

Figure 10 shows the results of uncertainty analysis for both the SD and PSO algorithm. The shaded region indicates the uncertainty in reservoir release during the period of this study.

The sources of uncertainty are related to input data measurement, model structure and parameter values. As shown in the figure, the width of the confidence interval in the PSO algorithm is significantly higher than the SD method. Therefore, the PSO algorithm has more uncertainty in the reservoir release values.

As mentioned before, the d-factor is used to determine the model uncertainties. Its values are 0.45 for the SD approach and 2.4 for the PSO algorithm. It can be said that the SD approach provides more appropriate planning in meeting the downstream demands compared to PSO.

Although the system dynamics method has many advantages in simulating complex systems, so far, the performance of the reservoir, regulatory dam, and diversion dam of the surface water control system has not been simulated simultaneously using this method. Therefore, in the present study, the performance of the system dynamics method in the management and simulation of the productivity of the water resource system for the desired integrated system was investigated. To evaluate the accuracy of this method, all results were compared with the PSO algorithm and the uncertainty of both methods was estimated. The following section provides an overview of the relevant results.

First, in order to validate the results of the SD model, the Dez Reservoir dam's performance was simulated for the statistical period of 2005–2014 using the above-mentioned model. The statistical indices of R², RMSE, and SE for the water volume behind the reservoir dam were 0.99, 1,690.3, and 0.18, respectively, and 0.98, 34.5 and, 0.08 for water released from the dam. Therefore, it can be said that the SD method has high accuracy in simulating the water volume of the Dez. Considering the capacity of the Sabili canal in the upstream of the regulatory dam, the amount of input water in this canal has been able to fully meet the demands of the canal. Also, the statistical indices of R², RMSE, and SE for released water from the regulatory dam were calculated as 0.976, 51.679, and 0.134, respectively.

The amount of water released from the regulatory dam was accumulated in the diversion dam, with the R², RMSE, and SE values of 0.85, 39.24, 0.18 for the main west canal and 0.86, 25.40, and 0.24 for the main east canal, and it has been able to meet the demands of the two mentioned canals to the desired extent. These results indicate the high accuracy of the SD model in the simultaneous simulation of the reservoir dam, regulatory dam, and diversion dam. In order to supply the demands of the Dez irrigation-drainage network, releases in the three cases of the real conditions, PSO algorithm optimization, and SD simulation were compared with one another based on which it can finally be said that, according to the obtained graphs as well as the deficiency values and the statistical indices, the release values simulated by the SD model were more in line with the downstream demands; and thus this model can have an effective role in the more accurate allocation of surface water resources in the region.

All relevant data are included in the paper or its Supplementary Information.