Optimal Parameter Estimation for Muskingum Model Using a CSS-PSO Method

Abstract

Limited availability of hydrologic data is a major hurdle for implementation of detailed hydrologic models. In cases where available data is limited, simple hydrologic models such as linear Muskingum model consisting of a minimum number (one or two) of model parameters are more desirable. As an alternative to the conventional mathematical approaches, this paper applies a new hybrid metaheuristic algorithm based on charged system search and particle swarm optimization for identifying the parameters of the linear Muskingum model. In order to evaluate the new algorithm, a numerical example is utilized and the results are compared to those of other algorithms. The results reveal the performance of the algorithm to optimize parameter estimation of the Muskingum model.

1. Introduction

Charged system search (CSS) is a relatively new metaheuristic optimization algorithm inspired by the governing laws of electrical physics and the Newtonian mechanics. In electrical physics, the electric charge can generate the electric field and exerts a force on other electrically charged objects. The electric field surrounding a point charge is specified by the laws of Coulomb and Gauss. Utilizing these principles, the CSS algorithm defines a number of solution candidates each of which is called charged particle (CP) and is treated as a charged sphere. Each CP can exert an electrical force on the other agents (CPs). These forces can change the position of other CPs according to Newton's second law. Finally, considering the Newtonian mechanics, the new positions of CPs are determined [1, 2]. Also a hybrid CSS and PSO (particle swarm optimization) algorithm is proposed by Kaveh and Talatahari [3]. The hybrid CSS-PSO has shown good performance in solving optimization problems in different areas such as optimal design of large-scale structures [4] and optimal design of laminated composite structures [5].

For the first time we apply this hybrid algorithm for optimal parameter estimation of Muskingum model. The Muskingum model is a popular model for flood routing which is basically an extension of the reservoir routing methods taking into account that the storage (S) in a channel is a function of both the inflow (I) and the outflow (O). This means that in addition to the prism storage there is a (positive or negative) wedge storage. However, the application of the model still suffers from an absence of an efficient method for parameter estimation. There are various mathematical techniques which have been used for estimating the parameters of the Muskingum model. These mathematical methods use gradient information to search the solution space near an initial starting point [6]. In order to evaluate the proposed method, a numerical example is utilized and the results are compared to those of other algorithms.

2. Charged System Search Algorithm

2.1. General Aspects

The charged system search (CSS) algorithm is a relatively novel metaheuristic based on the Coulomb and Gauss laws from electrical physics and the governing laws of motion from the Newtonian mechanics [1]. This algorithm can be considered as a multiagent approach, where each agent is a charged particle (CP). Each CP is considered as a charged sphere with radius a, having a uniform volume charge density and is equal to

q_{i} = \frac{fit (i) - fitworst}{fitbest - fitworst} i = 1,2, \dots, N,

(1)

where fitbest and fitworst are the best and the worst fitness of all the particles, fit (i) represents the fitness of the agent i, and N is the total number of CPs.

CPs can impose electric forces on the others, and their magnitude for the CP located in the inside of the sphere is proportional to the separation distance between the CPs, and for a CP located outside the sphere is inversely proportional to the square of the separation distance between the particles. The kind of the forces can be attractive or repelling and it is determined by using ar_ij, the kind of force parameter, as

a r_{i j} = {\begin{cases} + 1 & w . p . k_{t} \\ - 1 & w . p . 1 - k_{t}, \end{cases}

(2)

where ar_ij determines the type of the force, + 1 represents the attractive force and −1 denotes the repelling force, and k_t is a parameter to control the effect of the kind of force. Therefore, the resultant force is redefined as

\begin{matrix} F_{j} = q_{j} \sum_{i, i \neq j}^{} (\frac{q_{i}}{a^{3}} r_{i j} \cdot i_{1} + \frac{q_{i}}{r_{i j}^{2}} \cdot i_{2}) a r_{i j} p_{i j} (X_{i} - X_{j}) \\ j = 1,2, \dots, N, i_{1} = 1, i_{2} = 0 ⟺ r_{i j} < a \\ i_{1} = 0, i_{2} = 1 ⟺ r_{i j} \geq a, \end{matrix}

(3)

where F_j is the resultant force acting on the jth CP; r_ij is the separation distance between two charged particles defined as

r_{i j} = \frac{∥ X_{i} - X_{j} ∥}{∥ (X_{i} + X_{j}) / 2 - X_{best} ∥ + ∊},

(4)

where X_i and X_j are the positions of the ith and jth CPs, respectively, X_best is the position of the best current CP, and ∊ is a small positive number to avoid singularity. p_ij determines the probability of moving each CP toward the others as

p_{i j} = {\begin{cases} 1 & \frac{fit (i) - fitbest}{fit (j) - fit (i)} > rand \lor fit (j) < fit (i) \\ 0 & otherwise. \end{cases}

(5)

At the movement stage, each CP moves toward its new position under the action of the resultant forces and its previous velocity as

\begin{matrix} X_{j, new} = ran d_{j 1} \cdot k_{a} \cdot \frac{F_{j}}{m_{j}} \cdot Δ t^{2} + ran d_{j 2} \cdot k_{v} \cdot V_{j, old} \cdot Δ t + X_{j, old} \\ V_{j, new} = \frac{X_{j, new} - X_{j, old}}{Δ t}, \end{matrix}

(6)

where k_a is the acceleration coefficient; k_v is the velocity coefficient to control the influence of the previous velocity and here both are equal to 0.5; and rand_j1 and rand_j2 are two random numbers uniformly distributed in the range (0, 1). If each CP moves out of the search space, its position is corrected using the harmony search-based handling approach as described in [7]. In addition, to save the best design, a memory (charged memory) is utilized.

2.2. Hybrid PSO and CSS

The particle swarm optimization (PSO) involves a number of particles, which are initialized randomly in the space of the design variables. These particles fly through the search space and their positions are updated based on the best positions of individual particles and the best position among all particles in the search space which corresponds to a particle with the smallest weight. Also considering the position of a particle selected randomly from the swarm (passive congregation) can improve the performance of the algorithm.

The update moves a particle by adding a change velocity V_j^{k + 1} to the current position X_k^j as follows:

\begin{matrix} X_{j}^{k + 1} = X_{j}^{k} + V_{j}^{k}, \\ V_{j}^{k + 1} = ω V_{j}^{k} + c_{1} r_{1} (P_{j}^{k} - X_{j}^{k}) + c_{2} r_{2} (P_{g}^{k} - X_{j}^{k}) + c_{3} r_{3} (R_{j}^{k} - X_{j}^{k}), \end{matrix}

(7)

where ω is an inertia weight to control the influence of the previous velocity; r₁, r₂, and r₃ are three random numbers uniformly distributed in the range of (0, 1); c₁ and c₂ are two acceleration constants; c₃ is the passive congregation coefficient; P_j^k is the best position of the jth particle up to iteration k; P_g^k is the best position among all particles in the swarm up to iteration k; and R_k is a particle selected randomly from the swarm. In order to increase the exploration ability, in this paper the velocity of the particles is defined as

V_{j}^{k + 1} = ω V_{j}^{k} + c_{1} r_{1} (P_{j}^{k} - X_{j}^{k}) + c_{2} r_{2} (P_{g}^{k} - X_{j}^{k}) + c_{3} r_{3} (R_{j}^{k} - X_{j}^{k}) + c_{4} r_{4} (R d_{j}^{k} - X_{j}^{k}),

(8)

where c₄ is the exploration coefficient; r₄ is a uniformly distributed random number in the range of (0, 1); and Rd^k is a vector generated randomly from the search domain.

Compared to PSO, the CSS algorithm can be considered as a generalized form of the PSO containing its superiorities, but avoiding its disadvantages. Both are population-based algorithms and find optimum solutions by changing the position of the agents. However, the movement strategies are different for the CSS and PSO. The PSO algorithm utilizes the local best and the global best to determine the direction of the movement, while the CSS approach uses the governing laws from electrical physics and the Newtonian mechanics to determine the amount and the direction of a charged particle's movement. CSS can distinguish the local search and the global search phases and utilizes suitable relationships in these phases resulting in a good balance between the exploration and exploitation, while the greatest disadvantage of the PSO approach is the existence of some difficulties in controlling the balance between the exploration and exploitation.

In the present hybrid algorithm [3], the advantage of the PSO consisting of utilizing the local best and the global best is added to the CSS algorithm. The CM updating process is defined as follows:

C M_{i, new} = {\begin{cases} C M_{i, old} & fit (X_{i, new}) \geq fit (C M_{i, old}) \\ X_{i, new} & fit (X_{i, new}) \geq fit (C M_{i, old}) . \end{cases}

(9)

Considering the above new CM, the electric forces generated by agents are modified as

F_{j} = \sum_{i ∊ S_{1}}^{} (\frac{q_{i}}{a^{3}} r_{i j} \cdot i_{1} + \frac{q_{i}}{r_{i j}^{2}} \cdot i_{2}) a r_{i j} p_{i j} (C M_{i, old} - X_{j}) + \sum_{i ∊ S_{2}}^{} (\frac{q_{i}}{a^{3}} r_{i j} \cdot i_{1} + \frac{q_{i}}{r_{i j}^{2}} \cdot i_{2}) a r_{i j} p_{i j} (X_{i} - X_{j}),

(10)

where subtitles S₁ and S₂ denote two sets of the numbers which determine the number of the agents utilized to calculate the resultant force by employing the agents sorted in the CM and the current agents positions, respectively. If the coefficient k_i is defined as

k_{i} = (\frac{q_{i}}{a^{3}} r_{i j} \cdot i_{1} + \frac{q_{i}}{r_{i j}^{2}} \cdot i_{2}) a r_{i j} p_{i j},

(11)

then the resultant force formulae can be simplified as

F_{j} = k_{1} (C M_{g, old} - X_{j}) + k_{2} (C M_{j, old} - X_{j}) + \sum_{i ∊ S_{1}}^{} k_{i} (C M_{i, old} - X_{j}) + \sum_{i ∊ S_{2}}^{} k_{i} (X_{i} - X_{j}),

(12)

where the subtitle g denotes the number of the stored so far good position among all CPs. Therefore the first term directs the agents towards the global best position. When i = j, then the CM_{i, old} is treated similar to $P_{i}^{k}$ in the PSO as considered in the second term of the above equation. This will direct the agents towards the local best [4, 5].

3. Determining the Muskingum Model Parameters

3.1. Linear Muskingum Model

Hydrologic modeling plays a key role in water resources planning, development, and management. Limited availability of hydrologic data is a major hurdle for implementation of detailed hydrologic models. In cases where available data is limited, simple hydrologic models consisting of a minimum number (one or two) of model parameters are more desirable [8]. A simple conceptual model like the Muskingum model is very effective in the estimation of the hydrograph of the rise and fall of a river at any given point on the river during the course of a flood event.

The problem is solved by the techniques of flood routing, which is a process for calculating the shape of a flood wave along a river channel. McCarthy [9] analyzed the data from the Muskingum River in Ohio and found certain relationships: storage is proportional to outflow; and storage is also proportional to the difference between inflow and outflow [10]. In the linear Muskingum model, the following continuity and storage equations are used:

\begin{matrix} \frac{d S_{t}}{d t} \approx \frac{Δ S_{t}}{Δ t} = I_{t} - O_{t}, \\ S_{t} = K [x I_{t} - (1 - x) O_{t}], \end{matrix}

(13)

where S_t represents the storage volume within the reach at time t (L³); I_t and O_t are rates of the inflow and outflow at time t (L³ T⁻¹), respectively; K is the storage-time constant for the river reach (T), which has a value close to the flow travel time through the river reach; x is the dimensionless weighting factor commonly varying between 0.00 and 0.30 for stream channels and between 0.00 and 0.50 for reservoir storage; and Δt represents time step (T).

Choosing a suitable time interval for the routing period, the continuity equation in its finite difference form becomes

S_{i + 1} - S_{i} = \frac{(I_{i} + I_{i + 1}) Δ t}{2} - \frac{(O_{i} + O_{i + 1}) Δ t}{2}

(14)

using the linear Muskingum model:

O_{i + 1} = C_{1} I_{i} + C_{2} I_{i + 1} + C_{3} O_{i}, (i = 1,2, 3, \dots),

(15)

where

\begin{matrix} C_{1} = \frac{2 K x + Δ t}{Δ t + 2 K - 2 K x}, C_{2} = \frac{- 2 K x + Δ t}{Δ t + 2 K - 2 K x}, \\ C_{3} = \frac{2 K - 2 K x - Δ t}{Δ t + 2 K - 2 K x} \end{matrix}

(16)

in which O_i and I_i represent the outflow and inflow discharges at time t_i, respectively; and C₁, C₂, and C₃ are the Muskingum parameters where C₁ + C₂ + C₃ = 1. The routing procedure consists of the following steps.

Step 1. Assume the values of two parameters (K and x).

Step 2. Calculate the coefficients C₁, C₂, and C₃ using (16).

Step 3. Compute the magnitude of the outflow at the next time (O_{i + 1}) using (15).

Step 4. Repeat Step 3 for all times.

In the conventional method, the parameters K and x in the linear Muskingum model are graphically estimated by a trial and error procedure. After assuming x, the values of [XI_t + (1 – X)O_t] are computed using observed data in both upstream and downstream and then they are plotted against S_t. The value of x that minimizes the width of the plotted loop can be chosen as the correct value of x, and the slope of the straight line fitted through the loop gives K [11]. There are some disadvantages of the conventional trial and error method.

The satisfactory results obtained often need calculations of several times, rather than one time.

Because two parameters K, x are got through trying and error artificiality, the introduction of human errors will be unavoidable.

The results cannot be proved theoretically optimum solution of the Muskingum parameters and therefore cannot ensure the global optimum of C₁, C₂, and C₃.

This method not only increases the manual workload but also declines the accuracy of solution and it is difficult to programming.

Another traditional approach for determining Muskingum model parameters is the least square method (LSM). The LSM is a scientific solution for the Muskingum model under the condition of sufficient historical data. Because of its reliability and practicality, this method has been widely accepted. The biggest shortcoming of the LSM is the heavy dependence on historical data. If the historical data are complete and the time intervals are nearly equally distributed, the results of this approach will be good. On the contrary, if the historical data are incomplete, time intervals are too long or unequal, the results will be error greatly.

3.2. An Optimization Approach

In (15), C₁, C₂, and C₃ are hydrologic parameters of the model. The main challenge is to find a set of parameters that give a representative response. These parameters are determined using optimization techniques by obtaining a good match between the observed and calculated discharge. Choosing of a suitable objective function plays an important role in optimization problems. Here, we use the following objective function which is based on minimizing the sum of squares of deviations between observed and computed values of outflow:

SSQ = \sum_{j = 1}^{N_{o}} {(O_{o j} - O_{c j})}^{2},

(17)

where O_oj is observed value of the outflow at time t_j; O_cj is the computed value of the outflow at time t_j; j is an index varying from one to N_o, where N_o is the number of hydrograph ordinates of the observed hydrograph. Thus, the optimal parameter estimation for Muskingum model may be expressed as follows:

Find

{C_{1}, C_{2}, C_{3}}

(18)

To minimize

SSQ (C_{1}, C_{2}, C_{3})

(19)

Subject to

C_{1} + C_{2} + C_{3} = 1 .

(20)

4. Example of Hydrologic Parameter Estimation for Muskingum Model

To investigate the performance of applying the new algorithm to solve the parameter estimation problem of the linear Muskingum model, a typical problem is used as an example. The algorithms were implemented using MATLAB to run on a computer with an Intel Core i5 CPU, 2.53 GHz processor, and 3.00 GB RAM; also the penalty method is used to handle the constraint. In this paper, the model is applied to the south canal between Chenggou and Linqing rivers in China. The length of the river course (from Chenggou to Linqing) of South Canal in Haihe Basin reaches 83.8 km. There is no significant effect of the branches. The time interval in the calculation is taken as 12 hr. The data set from [13] is considered for illustration purpose. This example has been studied previously by Wang et al. [14, 15], Zahn and Xu [16], and Yang et al. [17] for testing different parameter estimation methodologies. Also Zhengxiang and Ling [12] applied and compared several intelligent algorithms on Muskingum routing model using this example. Therefore, the performance of the proposed algorithm can effectively be compared with the previously reported results obtained with this example.

In Tables 1 and 2, the optimum parameters and the objective function of the solutions obtained by the least squares method (LSM), GA, SA, PSO [12], ICA [6], and the proposed CSS-PSO are presented. For the proposed algorithm, a population of 10 particles is used. The results show that the SSQ is attained using new algorithm. It has been demonstrated that the new algorithm gets better results than those of other methods. The computed outflow hydrograph and the observed outflow hydrograph are shown in Figure 1. The plot depicts that the resulting outflow hydrograph, using the parameters estimated from the proposed method, closely follows the observed outflows.

Table 1:

Comparison of best SSQs from various methods.

	SSQ	C ₁	C ₂	C ₃

LSM [12]	1311.28	0.0843	0.4626	0.4531
GA [12]	1046.69	0.05476	0.45241	0.49283
SA [12]	1046.69	0.0547522	0.4524160	0.4928318
PSO [12]	1046.83	0.0547561	0.4524157	0.4928282
ICA [6]	1031.28	0.0100000	0.4774254	0.5116732
This study (CSS-PSO)	1020.22	0.0428869	0.4599482	0.4963921

Table 2:

Comparison of the observed and computed outflows (cm).

I (inflow)	Q_o (observed outflow)	Q_c (computed outflow)

261	228	228.00
389	300	303.29
462	382	379.73
505	444	440.58
525	490	481.83
543	513	511.44
556	528	532.89
567	543	549.16
577	553	562.31
583	564	572.02
587	573	578.94
595	581	586.22
597	588	591.10
597	594	593.61
589	592	591.17
556	584	574.44
538	566	556.44
516	550	536.62
486	520	512.04
505	504	507.29
477	483	492.86
429	461	462.43
379	420	422.26
320	368	373.05
263	318	319.87
220	271	271.24
182	234	227.79
167	193	197.69
152	178	175.21

Figure 1:

Observed and routed hydrographs.

5. Conclusion

In this research a powerful algorithm is developed by combining the CSS and PSO. The new algorithm in the CSS formulation introduces the PSO concept of local best and global best. The main advantages of the PSO consisting of directing the agents toward the global best (obtained by the swarm) and the local best (obtained by the agent itself) are added to the CSS algorithm to improve its performance. In the present approach, similar to the original CSS, each agent is affected by other agents considering the governing laws of electrical physics. However, the kind of the forces can be repulsive and attractive. For the proposed algorithm, when the calculations of the amount of forces are completed for the CPs and the new locations of agents are determined, the CM updating is performed.

Here, the hybrid CSS-PSO method has been proposed in order to improve the computational accuracy for parameter estimation of Muskingum model. The Muskingum method ignores the momentum equation and is based solely on the continuity equation. The proposed optimization approach found the best value of parameters in terms of the minimal sum of the square deviation between the observed and routed outflows. The performance of this approach was compared with those of other common methods such as GA, SA, PSO, and ICA. The results demonstrate that the new algorithm can achieve a high degree of accuracy to estimate the parameters. This method appears to offer good applicability in the hydrology field and further applications should be explored.

References

Kaveh

and Talatahari

, “A novel heuristic optimization method: charged system search,” Acta Mechanica, vol. 213, no. 3–4, pp. 267–289, 2010.

Kaveh

and Talatahari

, “Optimal design of skeletal structures via the charged system search algorithm,” Structural and Multidisciplinary Optimization, vol. 41, no. 6, pp. 893–911, 2010.

Kaveh

and Talatahari

, “Hybrid charged system search and particle swarm optimization for engineering design problems,” Engineering Computations, vol. 28, no. 4, pp. 423–440, 2011.

Kaveh

and Talatahari

, “Optimization of large-scale truss structures using modified charged system search,” International Journal of Optimization in Civil Engineering, vol. 1, no. 1, pp. 15–28, 2011.

Kaveh

Soudmand

B. H.

, and Sheikholeslami

, “Optimal design of laminated composite structures via hybrid charged system search and particle swarm optimization,” Asian Journal of Civil Engineering, vol. 14, no. 4, pp. 587–604, 2013.

Tahershamsi

and Sheikholeslami

, “Optimization to identify Muskingum model parameters using imperialist competitive algorithm,” International Journal of Optimization in Civil Engineering, vol. 1, no. 3, pp. 473–482, 2011.

Kaveh

and Talatahari

, “A charged system search with a fly to boundary method for discrete optimum design of truss structures,” Asian Journal of Civil Engineering, vol. 11, no. 3, pp. 277–293, 2010.

Ahmad

M. M.

Ghumman

A. R.

Ahmad

, and Hashmi

H. N.

, “Estimation of a unique pair of nash model parameters: an optimization approach,” Water Resources Management, vol. 24, no. 12, pp. 2971–2989, 2010.

McCarthy

G. T.

, “The unit hydrograph and flood routing,” in Proceedings of the North Atlantic Division Conference, U.S. Army Corps of Engineers, New London, Conn, USA, 1938.

10.

Geem

Z. W.

, “Parameter estimation of the nonlinear Muskingum model using parameter-setting-free harmony search,” Journal of Hydrologic Engineering, vol. 16, no. 8, pp. 684–688, 2011.

11.

Barati

, “Parameter estimation of nonlinear Muskingum models using nelder-mead simplex algorithm,” Journal of Hydrologic Engineering, vol. 16, no. 11, pp. 946–954, 2011.

12.

Zhengxiang

and Ling

, “Application and comparison of several intelligent algorithms on Muskingum routing model,” in Proceedings of the 2nd IEEE International Conference on Information and Financial Engineering (ICIFE '10), pp. 910–914, Chongqing, China, September 2010.

13.

G. J.

, “Parameter estimation for Muskingum models,” Journal of Hydrology, vol. 36, no. 3, pp. 40–43, 1977.

14.

Wang

Qiu

, and Xu

, “Hybrid chaotic genetic algorithms for optimal parameter estimation of Muskingum flood routing model,” in Proceedings of the IEEE International Joint Conference on Computational Sciences and Optimization (CSO '09), pp. 215–218, Hainan, China, April 2009.

15.

Wang

Kang

, and Qiu

, “Optimal parameter estimation for Muskingum model using a modified particle swarm algorithm,” in Proceedings of the 3rd IEEE International Joint Conference on Computational Sciences and Optimization (CSO '10), pp. 153–156, Anhui, China, May 2010.

16.

Zhan

and Xu

, “Application of ant colony algorithm to parameter estimation of Muskingum routing model,” Journal of Natural Disasters, vol. 14, no. 5, pp. 20–24, 2005.

17.

Yang

Jin

Chen

, and Wei

, “A new approach for parameter estimation of Muskingum routing model,” Journal of Catastrophology, vol. 13, pp. 1–6, 1998.