Modeling and Dynamical Analysis of the Water Resources Supply-Demand System: A Case Study in Haihe River Basin

Abstract

The relationship between water resources supply and demand is very complex and exhibits nonlinear characteristics, which leads to fewer models that can adequately manage the dynamic evolution process of the water resources supply-demand system. In this paper, we propose a new four-dimensional dynamical model to simulate the internal dynamic evolution process and predict future trends of water supply and demand. At the beginning, a new four-dimensional dynamical model with uncertain parameters is established. Then, the gray code hybrid accelerating genetic algorithm (GHAGA) is adopted to identify the unknown parameters of the system based on the statistic data (1998–2009). Finally, the dynamical analysis of the system is further studied by Lyapunov-exponent, phase portraits, and Lyapunov exponent theory. Numerical simulation results demonstrate that the proposed water resources supply-demand system is in a steady state and is suitable for simulating the dynamical characteristics of a complex water supply and demand system. According to the trends of the water supply and demand of several nonlinear simulation cases, the corresponding measures can be proposed to improve the steady development of the water resources supply-demand system.

1. Introduction

Water is indispensable to human survival, agriculture, and industry. Recently, there is an increasing shortage of water in most river basins mainly due to human activity and climate change [1, 2]. With the population growth, degradation of water quality, and loss of potential sources of freshwater, water demand is becoming increasing greater and water supply has been increasing scarce in many developing countries [3]. The internal dynamic evolution process and the relationship between water resources supply and demand are very complex and exhibiting nonlinear characteristics [4]. Water management is rapidly becoming more and more difficult due to the complex water resources system. Senge [5] discussed two types of complexity: (1) detail and (2) dynamic. Detail complexity is associated with systems with many components. Dynamic complexity, however, is in contact with effects that are separated by time or space. It is the dynamic complexity found in the water resources system that presents great difficulties for water management. Thus, how to analyse the internal dynamic evolved mechanism of the water resources supply-demand system and solve the issue of imbalanced water resources supply-demand has an important practical significance to the achievement of sustainable development and management of water resources.

The problem of water resources supply and demand has attracted much attention and a considerable amount of researches. Rayan [6] studied the water supply and demand of Sinai, Egypt, and introduced desalination options. Some previous works applied artificial neural networks (ANNs) to predict and forecast water resource system [7–9]. George et al. [10] proposed a city water balance model to analyse alternative strategies to manage the supply of water in order to come to terms with the water scarcity problems facing policy makers in Hyderabad. Qi and Chang [4] used system dynamics model to estimate domestic water demand under uncertain economic impacts from 2003 to 2009 for Manatee Country. Based on the data for meteorology, hydrology, soil, planting, vegetation, and socioeconomic development of the irrigation region in the middle reaches of the Heihe River basin, Ji et al. [11] established the model of balance of water supply and demand and assessed the security of water resources development. A probabilistic forecast model of future water demand for the city of Mecca in Saudi Arabia was proposed in [12]. Liu et al. [13] adopted WDF-ANN model to forecast water demand in Weinan city. Yurdusev et al. [14] predicted monthly water consumption by using feed-forward and radial-basis neural networks. However, the nonlinear dynamic behavior of the water resources supply-demand system is seldom considered by these methods.

Any very complicated dynamical system can be solved by a simple mathematical model [15]. The dynamical system method is one approach to nonlinear modeling. It has been applied increasingly in water environments and water resources system, such as hydrological time series prediction [16–18] and the dynamical analysis of runoff flows [19, 20]. In this paper, according to the current situation of the water resources supply and demand, we establish a new dynamical model to express the complex relationship between water resources supply and demand by using dynamical system theory. This model can be used to analyse the internal dynamic evolved mechanism of the water resources supply-demand system of a river basin or a city.

The outline of this paper is organized as follows: Section 2 introduces the study area and data. In Section 3, the water resources supply demand model is established and the method of parameter identification of the proposed model is introduced. Section 4 provides the results and discussions on the real water resources supply-demand system in Haihe River Basin. The conclusions of the paper are summarized in Section 5.

2. Study Area and Data

The HRB is the political, economic, and cultural center of China. However, water resources shortage in the HRB has become a serious problem in recent years due to the rapid economic development and climate change [21]. The location of the HRB is shown in Figure 1. The basin covers an area of 318200 km², with 60% being mountains in the western and northern part and 40% being plains in the eastern and southern parts. It belongs to continental monsoon climate zone and locates in a semihumid and semiarid region. The average annual precipitation is 524.69 mm, and the total amount of water resources is 374.39×10⁸m³ (1956–2008), the water volume per capita is 272 m³, merely 1/7 of the national average and 1/24 of the world average [22].

Figure 1 

The location of the Haihe River Basin.

The water resources demand and supply data from 1998 to 2009 are obtained from the water resources bulletin of HRB. Table 1 is the water resources statistic data (1998–2009) of HRB, where 𝑋 is water demand, 𝑌 is surface water supply, 𝑍 is groundwater supply, and 𝑊 is other kinds of water supply.

Table 1 

Water resources statistic data in Haihe River Basin (Unit: billion m³).

3. Modeling and Methods

3.1. Establishment of the Water Resources Supply-Demand Model

Let 𝑋(𝑡) be the demand for water resources, 𝑌(𝑡) the surface water supply (including diversion water supply), 𝑍(𝑡) the groundwater supply, and 𝑊(𝑡) other kinds of water supply (mainly including wastewater reuse, direct seawater utilization, and collection of rainwater). The model for the water resource system is established as follows: 𝑑𝑋𝑑𝑡=𝑎₁𝑋𝑋1−𝑀−𝑎₂𝑌−𝑎₃𝑍−𝑎₄𝑊,𝑑𝑌𝑑𝑡=−𝑏₁𝑌−𝑏₂𝑍+𝑏₃𝑋[𝑁−(𝑋−𝑍−𝑊)],𝑑𝑍𝑑𝑡=𝑐₁𝑋−𝑐₂𝑌−𝑐₃𝑍−𝑐₄𝑊,𝑑𝑊𝑑𝑡=𝑑₁𝑊𝑋−𝑑₂−𝑑₃𝑌, (1) where 𝑎_𝑖,𝑏_𝑖,𝑐_𝑖,𝑑_𝑖,𝑀,𝑁,𝐾 are constants, and 𝑀>0, 𝑁>0.

Additionally, 𝑎₁ is the elastic coefficient of the water resources demand; 𝑎₂, 𝑎₃, and 𝑎₄ are, respectively, the influence coefficients of the surface water supply, the groundwater supply, and other kinds of water supply on the water resources demand; 𝑏₁ is the influence coefficient of the surface water supply on the change rate of the surface water supply; 𝑏₂ is the influence coefficient of the groundwater supply on the surface water supply; 𝑏₃ is the influence coefficient of water resources demand on the surface water supply; 𝑐₁,𝑐₂,𝑐₃, and 𝑐₄ are, respectively, the influence coefficients of the water resources demand, the surface water supply, the groundwater supply, and other kinds of water supply on the groundwater supply; 𝑑₁ is the velocity constant of the recycled water supply; 𝑑₂ is the cost of recycled water supply; 𝑑₃ is the influence coefficient of the surface water supply on the other kinds of water supply; 𝑀 is the maximum value of water resources demand; 𝑁 is the valve value of diversion water supply.

In system (1), 𝑎₁𝑋(1−(𝑋/𝑀)) denotes that the rate of water resources demand 𝑑𝑋/𝑑𝑡 is in direct ratio to water resources demand when the current water resources demand of the study area is less than 𝑀, and 𝑑𝑋/𝑑𝑡 is inversely proportional to water resources demand when the water resources demand of the study area is larger than 𝑀; −𝑎₂𝑌, −𝑎₃𝑍, and −𝑎₄𝑊 denote that the increase in surface water supply, groundwater supply, and other kinds of water supply will reduce the demand of the water resources, respectively; −𝑏₁𝑌−𝑏₂𝑍 denotes that the rate of surface water supply will decrease with the surface water supply and groundwater supply; 𝑏₃𝑋[𝑁−(𝑋−𝑍−𝑊)] denotes that when the diversion water supply is less than the threshold value (i.e., 𝑋(𝑡)−𝑍(𝑡)−𝑊(𝑡)<𝑁), the rate of diversion water supply will increase with the increase of 𝑋(𝑡), and when diversion water supply is large enough (or 𝑋(𝑡)−𝑍(𝑡)−𝑊(𝑡)>𝑁), with economic development, the people and enterprises around the Yellow River and Yangtze River Basin will need a lot of water, which results in the rate of diversion water supply decreasing with the increase rate of 𝑋(𝑡); 𝑐₁𝑋 denotes that the increase in water resources demand will increase the groundwater supply; −𝑐₂𝑌−𝑐₃𝑍−𝑐₄𝑊 denotes that the decrease in surface water supply, the groundwater supply, and other kinds of water supply will decrease the rate of groundwater supply; 𝑑₁𝑊(𝑋−𝑑₂) denotes that the rate of reused water is in direct ratio to the amount of recycled water, which increases with the increase of water resources demand and decrease with the increase of the cost for reused water.

System (1) is an uncertainty system with unknown parameters, and it has different characteristics with different parameters. When a set of parameters are set as follows: 𝑎₁=0.06, 𝑎₂=0.2, 𝑎₃=0.09, 𝑎₄=0.15, 𝑏₁=0.22, 𝑏₂=0.2, 𝑏₃=0.2, 𝑐₁=0.13, 𝑐₂=0.14, 𝑐₃=0.12, 𝑐₄=0.2, 𝑑₁=0.16, 𝑑₂=0.2, 𝑑₃=0.001, 𝑀=1.1, and 𝑁=0.6, and the initial condition is set as [0.6,0.3,0.2,0.1], the Lyapunov exponents of the system (1) are 𝐿₁=−0.0156, 𝐿₂=−0.0371, 𝐿₃=−0.0853, and 𝐿₄=−0.1801. Because these Lyapunov exponents are all less than 0, the orbit of the system (1) is a limit cycle as shown in Figure 2. When 𝑎₁=0.06, 𝑎₂=0.08, 𝑎₃=0.3, 𝑎₄=0.2, 𝑏₁=0.05, 𝑏₂=0.04, 𝑏₃=0.04, 𝑐₁=0.1, 𝑐₂=0.2, 𝑐₃=0.05, 𝑐₄=0.04, 𝑑₁=0.2, 𝑑₂=0.3, 𝑑₃=0.001, 𝑀=1.1, and 𝑁=0.5, and the initial conditions are set as [0.7,0.2,0.1,0.1], the Lyapunov exponents of the system (1) are 𝐿₁=−0.0020, 𝐿₂=−0.0134, 𝐿₃=−0.0169, and 𝐿₄=−0.0592. The system has a stable focus point and its orbit is shown in Figure 3. When 𝑎₁=0.08, 𝑎₂=0.13, 𝑎₃=0.06, 𝑎₄=0.08, 𝑏₁=0.15, 𝑏₂=0.15, 𝑏₃=0.14, 𝑐₁=0.13, 𝑐₂=0.15, 𝑐₃=0.19, 𝑐₄=0.18, 𝑑₁=0.16, 𝑑₂=0.12, 𝑑₃=0.002, 𝑀=2.5, and𝑁=0.8, and initial condition is set as [0.8,0.1,0.3,0.2], the Lyapunov exponents of the system (1) are 𝐿₁=0.0421, 𝐿₂=−0.0104, 𝐿₃=−0.1144, and 𝐿₄=−0.2631. Since there are two Lyapunov exponents that are larger than zero, therefore, the system is hyperchaotic in this situation. System (1) has a chaotic attractor as shown in Figure 4.

(a)

(b)

(c)

Figure 2 

A limit cycle of system (1): (a) a 3D view (𝑋-𝑌-𝑍), (b) a 3D view (𝑋-𝑌-𝑊), and (c) the time series of 𝑋(𝑡),𝑌(𝑡),𝑍(𝑡),𝑊(𝑡).

(a)

(b)

(c)

Figure 3 

A stable focus point of system (1): (a) a 3D view (𝑋-𝑌-𝑍), (b) a 3D view (𝑋-𝑍-𝑊), and (c) the time series of 𝑋(𝑡),𝑌(𝑡),𝑍(𝑡),𝑊(𝑡).

(a)

(b)

(c)

Figure 4 

A chaotic attractor of system (1): (a) a 3D view (𝑋-𝑌-𝑍), (b) a 3D view (𝑋-𝑌-𝑊), and (c) the time series of 𝑋(𝑡),𝑌(𝑡),𝑍(𝑡),𝑊(𝑡).

3.2. Parameter Identification

To illustrate the behavior of water resources supply-demand system in a specified area, the unknown parameters should be identified first. Parameter identification can be seen as a global optimization problem. Many methods for the global optimization problem have been proposed, which can be roughly divided into two categories: deterministic optimization methods and stochastic optimization methods. The deterministic approaches take advantage of the analytical properties of the problem to generate a sequence of points that converge to a global optimal solution and sometimes can give a guaranteed solution [23, 24]. However, the deterministic optimization may be trapped in the local extreme point when the dimension is high and there are numerous local optima. Sometimes, traditional deterministic methods cannot obtain the global optimization efficiently due to its complex computation and strict application conditions [25]. Additionally, other research results suggest that the deterministic global algorithms have lower efficiency for the optimization problems with a large number of parameters [26]. Therefore, the traditional deterministic global algorithms are not considered to solve our problem which not only contains many parameters, but also shows high nonlinearity and complex dynamic characteristics.

Genetic algorithm (GA) is a type of stochastic optimization method based on the mechanics of natural selection and natural genetics. It can provide a robust procedure not only to explore broad and promising regions of solutions, but also to avoid being trapped in the local optimization [27]. Many researchers have applied the GA techniques to identify linear and nonlinear systems. Based on many improved genetic algorithms, Yang et al. [28] proposed a gray-encoded hybrid accelerating genetic algorithm (GHAGA) with Nelder-Mead simplex searching operator and simplex algorithm for the global optimization of dynamical systems. The GHAGA method is an improved genetic algorithm. Compared with traditional genetic algorithm and pure random search algorithm, the GHAGA method has the advantage of higher calculation precision and efficiency, rapider convergent speed, and wider scope of application. It has been proved as an effective method to be used in various complex water environment optimization problems [28, 29]. Therefore, the GHAGA method is chosen to identify the parameters of the water resources supply-demand system (1).

3.2.1. Discretization

Before determining the parameters, system (1) needs to be discretized. Using the discrete transform, the water resources system (1) is discretized into system (2): 𝑋()(𝑘)𝑎𝑘+1=𝑋+𝑡₁𝑋(𝑘)𝑋(𝑘)1−𝑀−𝑎₂𝑌(𝑘)−𝑎₃𝑍(𝑘)−𝑎₄𝑊(𝑘,𝑌()(𝑘))𝑘+1=𝑌+𝑡−𝑏₁𝑌(𝑘)−𝑏₂𝑍(𝑘)+𝑏₃𝑋(𝑘(𝑋(𝑘)(𝑘)(𝑘,𝑍()(𝑘)𝑐)[𝑁−−𝑍−𝑊))]𝑘+1=𝑍+𝑡₁𝑋(𝑘)−𝑐₂𝑌(𝑘)−𝑐₃𝑍(𝑘)−𝑐₄𝑊(𝑘,𝑊()(𝑘)𝑑)𝑘+1=𝑊+𝑡₁𝑊(𝑘)𝑋(𝑘)−𝑑₂−𝑑₃𝑌(𝑘). (2)

3.2.2. The Algorithm of GHAGA

The GHAGA with Nelder-Mead simplex searching operator and simplex algorithm is used for the parameter identification. A flowchart for the identification of the system parameters by GHAGA is given in Figure 5, and the detailed steps of the algorithm are given in the following.

Figure 5 

The flowchart for the identification of the system parameters by GHAGA.

Consider the parameter optimization of the water resources model as min𝑓𝑥₁,𝑥₂,…,𝑥_𝑛s.t.𝑎_𝑖≤𝑥_𝑖≤𝑏_𝑖,for𝑖=1,2,…,𝑛, (3) where 𝑥={𝑥_𝑖,𝑖=1,2,…,𝑛}, 𝑥_𝑖 is a parameter to be optimized, 𝑓 is an objective function, and 𝑓≥0, for system (1), 𝑛=16.

Step 1 (gray encoding). Suppose gray encoding length is 𝑒 in every parameter, the 𝑖th parameter range is the interval [𝑎_𝑖,𝑏_𝑖], and then each interval is divided into 2^𝑒−1 subintervals: 𝑥_𝑖=𝑎_𝑖+𝐼_𝑖⋅𝑐_𝑖, (4) where the length of subinterval of the 𝑖th parameter 𝑐_𝑖=(𝑏_𝑖−𝑎_𝑖)/(2^𝑒−1) is constant. The searching location 𝐼_𝑖 is a nonnegative integer, and 0≤𝐼_𝑖<2^𝑒 for 𝑖=1,2,…,𝑛. The gray code array of the 𝑖th parameter is denoted by the grid points of {𝑑(𝑖,𝑘)∣𝑘=1,2,…,𝑒} for every individual: 𝐼_𝑖=^𝑒_𝑗=1^𝑒_𝑘=𝑗𝑑(𝑖,𝑘)⋅2^𝑗−1, (5) where ⨁ denotes the operator of addition modulo 2 on {0,1} [29]. The GHAGA uses a code of parameters instead of the parameters and works on a population of points and not on one single point.

Step 2 (randomly generating the initial population). Initially, the chromosomes are generated at random in gray-encoded genetic algorithm, and p-chromosomes in father population are 𝐼_𝑖𝑗=int𝑢𝑖,𝑗⋅2^𝑒,𝑖=1,2,…,𝑛,𝑗=1,2,…,𝑝, (6) where 𝑢(𝑖,𝑗) is uniform random numbers, 𝑢(𝑖,𝑗)∈[0,1], 𝐼_𝑖(𝑗) is a searching location, and int(⋅) is an integer function. From (3), the p-corresponding chromosomes are 𝑑(𝑖,𝑘,𝑗), for 𝑖=1,2,…,𝑛; 𝑘=1,2,…,𝑒; 𝑗=1,2,…,𝑝. To cover homogeneously the whole solution space and to avoid the risk of having too much individuals in the same region, a large uniformity random population is selected in this algorithm.

Step 3 (decoding and fitness evaluation). Decoding of 𝑑(𝑖,𝑘,𝑗) for (𝑖=1,2,…,𝑛;𝑘=1,2,…,𝑒;𝑗=1,2,…,𝑝) works through (4) and (5) and then corresponding parameter 𝑥_𝑖(𝑗) is obtained. Substitution of 𝑥_𝑖(𝑗) into (1) produces the objective function 𝑓(𝑗). The smaller the value 𝑓(𝑗) is, the higher the fitness of its corresponding 𝑖th chromosome will be. Then, the fitness function of 𝑖th chromosome is defined: 𝐹𝑗=1𝑓𝑗².+0.1 (7)

Step 4 (selection). The chromosomes in the initial father population are selected by a known probability 𝑃_𝑠(𝑗)∶𝑃_𝑠(𝑗)=𝐹(𝑗)/∑^𝑛_𝑖=1𝐹(𝑖). Such two groups of p-chromosomes are selected by the above probabilities.

Step 5 (two-point crossover). Perform crossover on each chromosome pair according to probability 𝑃_𝑐 to generate two offsprings. For two-point crossover, two crossing points 𝐼₁=int(𝑈₁⋅(𝑒+1)) and 𝐼₂=int(𝑈₂⋅(𝑒+1)) are randomly chosen, where 𝑈₁,𝑈₂ are uniformity random numbers, where {𝑈₁,𝑈₂}∈[0,1]. The two-point crossover between two individuals is performed through the crossing probability 𝑃_𝑐: 𝑑^′₁=𝑑𝑖,𝑘,𝑗²𝑖,𝑘,𝑗,𝑘∈𝐼₁,𝐼₂,𝑑₁𝑖,𝑘,𝑗,𝑘∉𝐼₁,𝐼₂;𝑑^′₂=𝑑𝑖,𝑘,𝑗¹𝑖,𝑘,𝑗,𝑘∈𝐼₁,𝐼₂,𝑑₂𝑖,𝑘,𝑗,𝑘∉𝐼₁,𝐼₂. (8)(9) In order to enhance the diversity of population, the crossing probability is set as 𝑃_𝑐≥0.5.

Step 6 (two-point mutation). The operator of two-point mutation is for four random numbers {𝑉₁,𝑉₂,𝑉₃,𝑉₄}∈[0,1]. If 𝑉₁≤0.5, the offspring is computed by (8). Otherwise, the offspring 𝑑^′(𝑖,𝑘,𝑗) is computed by (9). Let the mutating probability be 𝑝_𝑚∈[0,1], the two points be 𝐽₁=int(𝑉₂⋅(𝑒+1)), 𝐽₂=int(𝑉₃⋅(𝑒+1)), and then the above offspring is mutated from 𝑑₃𝑖,𝑘,𝑗_{(𝑖,𝑘,𝑗)}=0𝑘=𝐽₁,𝑉₄≤𝑝_𝑚,𝑑^′𝑖,𝑘,𝑗_{(𝑖,𝑘,𝑗)}=11𝑘=𝐽₁,𝑉₄≤𝑝_𝑚,𝑑^′𝑖,𝑘,𝑗_{(𝑖,𝑘,𝑗)}=00𝑘=𝐽₂,𝑉₄≤𝑝_𝑚,𝑑^′𝑖,𝑘,𝑗_{(𝑖,𝑘,𝑗)}=11𝑘=𝐽₂,𝑉₄≤𝑝_𝑚,𝑑^′𝑖,𝑘,𝑗_{(𝑖,𝑘,𝑗)}𝑑=0_′𝑖,𝑘,𝑗_{(𝑖,𝑘,𝑗)}otherwise, (10) where a new offspring 𝑑₃(𝑖,𝑘,𝑗) can be computed by a mutating probability 𝑝_𝑚.

Step 7 (simplex evolution and alternating). The Nelder-Mead algorithm is a useful, local descent algorithm, which does not make use of the objective function derivatives [30]. The best point in the previous phase becomes a new initial solution in the Nelder-Mead simplex algorithm, and then a new best point is obtained by this Nelder-Mead simplex algorithm. The new best point inside the offspring will be inserted to replace the worst one in the previous phase. Repeat Step 3 to Step 7 until the evolution times 𝑄 is met.

Step 8 (accelerating cycle). The parameter ranges of 𝑛_𝑠-excellent individuals obtained by 𝑄-times of the Nelder-Mead simplex evolution alternating become the new ranges of the parameters and then the whole process back to the Gray-encoding. The GHAGA computation is over until the algorithm running times get to the design 𝑇 times or there exists a chromosome 𝑐_fit whose fitness satisfies a given criterion. In the former case the 𝑐_fit is the fittest chromosome or the most excellent chromosome in the population. The chromosome 𝑐_fit represents the solution.

4. Results and Discussions

4.1. The Results of Parameters Identification

The GHAGA is calculated by using Matlab. Table 2 presents the variable initial intervals of the 16 unknown parameters in system (1) which are given according to the meaning of the parameters and the empirical research. The parameters of the GHAGA are selected as follows: the length 𝑒=10, the population size 𝑛=300, the number of excellent individuals 𝑛_𝑠=10, the times of simplex evolution alternating 𝑄=5, the crossover probability 𝑝_𝑐=1.0, the mutation probability 𝑝_𝑚=0.5, and the times of simplex searching 𝑚=600. The computational results of parameter identification for system (1) are also listed in Table 2.

Table 2 

The parameters’ initial intervals and identified values. P stands for parameters, II stands for initial intervals, and IP stands for identified parameters.

4.2. The Dynamical Analysis of the Real System

When the parameters values are taken as shown Table 2, system (1) has four equilibrium points: 𝑂(0,0,0,0), 𝑆₁(0.9093,0.7669,0.8118,0.3969), 𝑆₂(0.8856,0.4171,3.6269,−0.9523), and 𝑆₃(0.3586,−2.8007,6.1825,0.0527). By calculation, the eigenvalues of the Jacobi matrix of the system (1) at 𝑂(0,0,0,0) are 𝜆₁=0.6921>0, 𝜆₂=0.0320>0, 𝜆_2,3=−0.1155>0, and 𝜆₄=−0.0876>0. Because four eigenvalues are not all negative, therefore 𝑂(0,0,0,0) is an unstable point and the system is unstable at the point 𝑂(0,0,0,0).

The eigenvalues of the Jacobi matrix of the system (1) at 𝑆₁(0.9093,0.7669,0.8118,0.3969) are 𝜆₁=−0.0022<0, 𝜆₂=−0.5721<0, and 𝜆_3,4=−0.0519±0.0677𝑖. Obviously, the real parts of all the four eigenvalues are negative, so the water resources supply-demand system is stable when the demand of water resources in Haihe River Basin, the surface water supply, the groundwater supply, and the other kinds of water supply are in the local area of the equilibrium point 𝑆₁(0.9093,0.7669,0.8118,0.3969). It indicates that water resources supply-demand system has the capacity of self-adjustment and demonstrates the status of steady development in the local area of 𝑆₁(0.9093,0.7669,0.8118,0.3969).

Similarly, we can prove that the water resources supply-demand system is unstable when the demand of water resources in Haihe River Basin, the surface water supply, the groundwater supply, and the other kinds of water supply are in the local area of the equilibrium points 𝑆₂ and 𝑆₃.

4.3. Numerical Simulations

In this section, based on the results of parameter identification by GHAGA, the numerical simulation results for the real water resources supply-demand system are obtained, and corresponding measures are given to ensure the steady development of the water resources supply-demand system.

Case 1. The result of parameter identification is used as the system parameters (Table 2) and the initial condition is set as [0.7,0.3,0.1,0.5], and the water resources supply-demand system is stable as shown in Figure 6. Under the current initial situation and reasonable water resources demand, the relationship between water resources’ supply and demand of the HRB will be in a stable state. From Figure 6, we can see that the trend of the water resources demand will increase firstly and then decrease; the surface water supply will be in a growth trend until steady state; the groundwater supply will increase firstly and then decrease until a steady state; other kinds of water supply will be in a growth state. According to the above analysis of the trend of the three kinds of water supply sources, we can obtain that, in order to maintain the balance between water resources supply and demand, we need to save water and improve water use efficiency, give priority to use of south-to-north water, and improve the use of recycled water.

(a)

(b)

(c)

(d)

Figure 6 

The time series of 𝑋(𝑡),𝑌(𝑡),𝑍(𝑡),𝑊(𝑡).

Case 2. Maintain other parameters constant and initial conditions are set as [0.7,0.3,0.1,0.5], when the maximum value 𝑀>1.5, the time series 𝑋(𝑡), 𝑌(𝑡), 𝑍(𝑡), and 𝑊(𝑡) are shown in Figure 7. From Figure 7, we can obtain that the values of 𝑌(𝑡), 𝑍(𝑡), and 𝑊(𝑡) are too large and beyond the scope of practical significance. It indicates that if the maximum demand of water resources is set to be very large; that is to say, if the water resources demand is not controlled, the system will be in an unstable state. Therefore, we must save water and control water demand.

(a)

(b)

(c)

(d)

Figure 7 

The time series of 𝑋(𝑡),𝑌(𝑡),𝑍(𝑡),𝑊(𝑡) when 𝑀>1.5.

Case 3. Set 𝑑₂=1.4 and maintain other parameters constant and the initial conditions are set as [0.9,0.5,0.7,0.3], the simulation result is shown in Figure 8. The figure shows that the trend of the other kinds of water resources supply will decrease until zero; the groundwater supply will increase until larger than 1, and this will lead to overexploitation of the groundwater supply. Therefore, with the cost of recycled water supply growing, the other kinds of water supply will decrease. That is to say, other types of water supply are dependent on the cost of the recycled water supply.

(a)

(b)

(c)

(d)

Figure 8 

The time series of 𝑋(𝑡),𝑌(𝑡),𝑍(𝑡),𝑊(𝑡) when 𝑑2=1.4.

Case 4. Maintain other parameters constant and vary valve value 𝑁, the water resources supply-demand system still maintains a stable state. It indicates that the valve value 𝑁 has no influence on the system. However, if the valve value 𝑁 is set to be very large, it will cause water resources shortage in other regions and further affect the economic development of the other regions. Therefore, an appropriate valve value needs to be selected.

5. Conclusions

In this paper, a novel model is proposed to analyse the relationship between water resources supply and demand. The main conclusions are as follows.

(1) Using dynamical system theory, we establish a new four-dimensional water resources supply-demand system with unknown parameters. This model can be used to analyse the complex nonlinear characteristics among the water resources demand, surface water supply, groundwater supply, and other kinds of water supply of a river basin or a city.

(2) Selecting the Haihe River Basin as a study area, based on its statistic data (1998–2009), we identify the unknown parameters of the system with the gray code hybrid accelerating genetic algorithm (GHAGA). Therefore, a real water resources supply-demand system of the Haihe River Basin is determined.

(3) Numerical simulation results provide the trend of the water resources demand and supply. The results show that the obtained water resources supply-demand system has the capacity of self-regulatory and demonstrates a steady state. Several cases are analysed and contribute to the formulation of policy and the management of water resources. Finally, the corresponding measures are proposed to ensure the steady development of the water resources supply-demand system.

Our research provides a new method for water resources system analysis based on dynamical system theory. This method can also be used to discuss the balance between water resources supply and demand in other regions or cities.

Conflict of Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.

Acknowledgments

This work was supported by the National Basic Research Program of China (no. 2010CB951104), the the Funds for Creative Research Groups of China (no. 51121003), and the Project of National Natural Foundation of China (no. 51379013).