Calibrating RZWQM2 model using quantum-behaved particle swarm optimization algorithm

Computers and Electronics in Agriculture 113 (2015) 72–80

Contents lists available at ScienceDirect

Computers and Electronics in Agriculture journal homepage: www.elsevier.com/locate/compag

Calibrating RZWQM2 model using quantum-behaved particle swarm optimization algorithm Maolong Xi a,⇑, Zhiming Qi b, Ye Zou a, G.S. Vijaya Raghavan b, Jun Sun c a

Department of Control Technology, Wuxi Institute of Technology, Wuxi, Jiangsu 214121, China Department of Bioresource Engineering, McGill University, Sainte-Anne-de-Bellevue, QC H9X 3V9, Canada c Department of Computer Science, Jiangnan University, Wuxi, Jiangsu 214122, China b

a r t i c l e

i n f o

Article history: Received 28 April 2014 Received in revised form 9 January 2015 Accepted 3 February 2015 Available online 19 February 2015 Keywords: QPSO Calibration RZWQM2 Single-objective Multiple-objective

a b s t r a c t The RZWQM2 model employs a number of parameters and coefficients requiring calibration and validation to accurately predict how interacting environmental conditions and management approaches influence crop development, water flow and pollutant levels in an agricultural production system. Given these requirements, an auto-calibration tool capable of optimizing RZWQM2 parameters would greatly enhance model efficiency over manual trial-and-error calibration. With a single parameter controlling convergence speed, the easy-to-implement quantum-behaved particle swarm optimization (QPSO) algorithm has performed well in deriving solutions to a wide range of optimization problems. Five years (2005–2009) of yield, drain flow, and NO 3 —N loss data from a subsurface-drained corn-soybean field in Iowa were employed in assessing the feasibility of using such an algorithm as an adjunct to RZWQM2. Using the QPSO algorithm, RZWQM2 parameters were calibrated using single- and multiple-objective functions. The algorithm-enhanced and standard models’ performances in simulating yield, drainage, and NO3–N loss were compared on the basis of the percent bias (PBIAS), Nash–Sutcliffe efficiency coefficient (NSE), and the ratio of root mean squared error to standard error (RSR). With the QPSO algorithm the model calibration by individual parameter yielded PBIAS, NSE and RSR values of ±10%, over 0.90, and less than 0.44, respectively, while for the collectively calibrated model values were ±12%, over 0.83, and under 0.41, comparable to the performance achieved through manual calibration. This study suggests that the models calibrated with the QPSO algorithm operated in a satisfactory manner, and this method can be used successfully in estimating parameters in RZWQM2. Ó 2015 Elsevier B.V. All rights reserved.

1. Introduction A typical swarm intelligence (SI) algorithm, particle swarm optimization (PSO) algorithm, originally developed by Kennedy and Eberhart (1995), was inspired by the social behavior of bird flocks and fish schools. Initialized with a population of random solutions, a PSO system searches for optima by updating generations through a heuristic, stochastic, parallel, and distributed optimization method based on swarm behavior (Clerc and Kennedy, 2002). Having no evolution operators (e.g., crossover and mutation) such as exist in genetic algorithms (GA), the PSO is independent of the problem’s internal structures. It therefore rarely requires the analysis of objective function properties, making it a good choice for optimization problems where the objective function is only implicitly given or is not differentiable. ⇑ Corresponding author. Tel.: +86 0510 81838727. E-mail addresses: [email protected], [email protected] (M. Xi). http://dx.doi.org/10.1016/j.compag.2015.02.002 0168-1699/Ó 2015 Elsevier B.V. All rights reserved.

Comparable in performance to evolutionary algorithms, the PSO algorithm may therefore be considered an alternative to such methods (Angeline, 1998; Poli et al., 2007; Kuok and Chan, 2012). Motivated by quantum mechanics and trajectory analysis of PSO (Clerc and Kennedy, 2002), Sun et al. (2004a) used a strategy based on a quantum-delta-potential-well model to sample around the previous best points. The mean best position was incorporated into the algorithm and the newly developed PSO termed a quantum-behaved particle swarm optimization (QPSO) model (Sun et al., 2004b). With a probabilistic algorithm shown to exhibit global convergence (Sun et al., 2012), the QPSO differs significantly from PSO in its iterative equation and in not requiring a velocity vector for particles. Having only one parameter to adjust, the QPSO’s implementation is easier than that of classical gradientbased techniques. QPSO algorithms have been implemented in the design of an infinitesimal dipole model for dielectric resonator antennas (Mikki and Kishk, 2006), in parameter optimization and model identification (Gao and Tong, 2006; Xi et al., 2007; Sabat

73

M. Xi et al. / Computers and Electronics in Agriculture 113 (2015) 72–80

et al., 2009), as well as applications in communication networks (Sun et al., 2006), control engineering (Xi et al., 2006), fuzzy systems (Tang and Xue, 2008), image processing (Lei and Fu, 2008), and energy management (Coelho and Mariani, 2008). An agricultural system model widely used in simulating plant growth, water movement, and the fate of nutrient and pesticides, the Root Zone Water Quality Model (RZWQM2), like most agricultural models, requires many of its input parameters to be calibrated in order to perform well (Ahuja et al., 2000). Many strategies, including automatic calibration methods, have been applied in an effort to raise the efficiency and accuracy of model calibration: e.g., a GA was used to determine cultivar coefficients for crop simulation models (Pabico et al., 1999), and simulated annealing was implemented in a soil water balance model (Calmon et al., 1999). Notwithstanding its time-consuming and somewhat subjective nature, and the fact that model performance depends largely on user’s experience, a traditional trial-and-error calibration method remains that most often employed in calibrating the RZWQM2 model (Thorp et al., 2007; Qi et al., 2011; Ma et al., 2012). However, a systematic and repeatable approach for RZWQM2 parameter estimation was achieved using PEST parameter estimation software (Nolan et al., 2010; Fang et al., 2010; Malone et al., 2010). The software implemented a Gauss– Marquardt–Levenberg nonlinear scheme, a gradient-based optimization strategy that combines the Gauss–Newton algorithm and a gradient descent method. The PEST method provides a numerical solution to the mathematical problem of minimizing a sum of weighted least squares between model outcomes and corresponding field data. It can produce a fast and efficient convergence to an objective function minimum. However, when objective functions are discontinuous or characterized by many local optima, the PEST search procedure can be sensitive to the choice of initial values. Furthermore, PEST method is integrated in software whose objective function format is fixed and cannot be modified freely by users to address different optimization problems. Therefore, other approaches to optimizing RZWQM2 parameters should be investigated. Since, to the best of our knowledge, the application of swarm intelligence algorithms for calibration of the RZWQM2 model has not been reported elsewhere, the goal of this work was to employ and assess the efficacy of a QPSO algorithm in calibrating RZWQM2.

c1, c2 are acceleration coefficients which control how far a particle will move in a single iteration; they are often assigned a value of 2, r1  U(0, 1), r2  U(0, 1) are elements from two uniform random sequences in the range of (0, 1), w is the inertial weight, akin to the momentum term in a gradient descent neural-network training algorithm; it is often set to decrease from 0.9 to 0.4 during the optimization process (Shi and Eberhart, 1998), Pi,j = (Pi,1, Pi,2,   , Pi,n) is the local best position of the particle, Pg,j = (Pg,1, Pg,2,   , Pg,n) is the global best position of the swarm, and Vi,j = (Vi,1, Vi,2,   , Vi,n) is the current velocity of the particle. In a previous study, a trajectory analysis under the PSO algorithm, Clerc and Kennedy (2002) demonstrated that each particle must converge to its local attractor pi,j to guarantee convergence of the algorithm. These coordinates are defined as:

pi;j ðtÞ ¼

u1 Pi;j ðtÞ þ u2 P g;j ðtÞ ; where u1 ; u2  Uð0; 1Þ u1 þ u2

ð3Þ

or

pi;j ðtÞ ¼ u  Pi;j ðtÞ þ ð1  uÞ  Pg;j ðtÞ; where u  Uð0; 1Þ

ð4Þ

Inspired by this analysis of the convergence of the traditional PSO and assuming a one-dimensional delta potential well where in each particle exhibits quantum behavior for each dimension at point pi, Sun et al. (2004a) solved the Schrödinger equation for each dimension. They thereby obtained the normalized probability density function Q and distribution function F for each component of the particle’s current position:

  Q X i;j ðt þ 1Þ ¼

1 2jpi;j ðtÞXi;j ðtþ1Þj=Li;j ðtÞ e Li;j ðtÞ

ð5Þ

  F X i;j ðt þ 1Þ ¼ e2jpi;j ðtÞXi;j ðtþ1Þj=Li;j ðtÞ

ð6Þ

where Li,j(t) is the standard deviation of the distribution, which determines the search scope of each particle. Employing the Monte Carlo method, the position of the particle can then be obtained as:

Li;j ðtÞ lnð1=uÞ 2

2. Materials and methods

X i;j ðt þ 1Þ ¼ pi;j ðtÞ 

2.1. Quantum-behaved particle swarm optimization

where u is a random number uniformly distributed in the interval (0, 1). A global point called the Mainstream Thought or Mean Best Position, m(t), of the population was introduced into PSO by Sun et al. (2004b). The global point is defined as the average of the pbest positions of all particles:

Having few parameters to be adjusted along with the ability to quickly and cheaply generate accurate solutions to optimization problems, PSO has gained in popularity. Rather than using evolutionary algorithms’ evolutionary operators to manipulate individuals, PSO relies on an exchange of information between individuals. Particles, representing problem solutions, fly over the problem landscape following their own path and that of the current best particles. Each particle is attracted to the best location that it has found, and also the best location found by other particles in its topological neighborhood within the problem space. The equations for assigning the particle’s velocity and for changing its positions are:

V i;j ðt þ 1Þ ¼ w  V i;j ðtÞ þ c1  r 1 ðtÞ  ½Pi;j ðtÞ  X i;j ðtÞ þ c2  r 2 ðtÞ    Pg;j ðtÞ  X i;j ðtÞ

ð1Þ

X i;j ðt þ 1Þ ¼ X i;j ðtÞ þ V i;j ðt þ 1Þ

ð2Þ

where

mðtÞ ¼ ðm1 ðtÞ; m2 ðtÞ;    ; mn ðtÞÞ ¼

M M M 1X 1X 1X Pi;1 ðtÞ; Pi;2 ðtÞ;    ; Pi;n ðtÞ M i¼1 M i¼1 M i¼1

ð7Þ

! ð8Þ

where M is the population size, and Pi is the pbest position of particle i. The values of Li,j(t) and the position Xi,j(t + 1) are evaluated as:

  Li;j ðtÞ ¼ 2b  mj ðtÞ  X i;j ðtÞ   X i;j ðt þ 1Þ ¼ pi;j ðtÞ  b  mj ðtÞ  X i;j ðtÞ  lnð1=uÞ

ð9Þ ð10Þ

74

M. Xi et al. / Computers and Electronics in Agriculture 113 (2015) 72–80

where b is the Contraction–Expansion coefficient, the only parameter that can be tuned to control the algorithms’ convergence speed. An algorithm which incorporates Eqs. (3), (8) and (10) is termed a QPSO. The QPSO algorithm’s procedure is shown in Fig. 1. 2.2. Root Zone Water Quality Model (RZWQM2) RZWQM2 (version 2.00.2010) is a one-dimensional agricultural system model consisting of components for hydrology, nutrition and pesticide transport and transformation, plant growth, and management practices (Ahuja et al., 2000; Ma et al., 2005, 2006). Infiltration from rainfall, irrigation, or snow melt is computed by a modified Green–Ampt approach, and the water redistribution in the soil profile is simulated by the Richards equation and considers surface evaporation and plant uptake as sinks. Water stored in the soil profile, if in excess of the field capacity, is drained to build up a water table above the impermeable layer and the subsequent tile drainage flux is calculated using the steady state Hooghoudt equation. In addition to subsurface drainage, lateral flow and deep seepage are quantified by user-defined parameters for constant bottom layer water flux rate and the hydraulic gradient, respectively. For the upper boundary, soil evaporation and plant transpiration are estimated by the double layer model of Shuttleworth and Wallace (1985), which is an extension of the Penman–Monteith concept. The nutrient chemistry processes model incorporated in RZWQM2 is OMNI (Shaffer et al., 2000), a state-of-the-art model for carbon and nitrogen cycling in soils. The coupled DSSAT family

of crop growth models (Jones et al., 2003) enhances the capability of RZWQM2 to describe crop establishment and water and nutrient uptake. The parameters for crop cultivars commonly planted in the United States have been tested and are well documented for the users of the model. Users can conveniently adjust the crop root water uptake by changing the soil root growth factor for each calculated soil layer. 2.3. Field data and parameters A field sub-data set collected from a subsurface drained corn– soybean field in Iowa was used to examine the feasibility of using QPSO in the optimization of RZWQM2 parameters. Measured data for crop grain yield, subsurface drainage flow, and NO 3 —N loss from four drainage plots, labeled 14-3, 24-2, 18-1, and 10-2, were extracted from Qi et al. (2011). The experiment was conducted in 2005–2009 in fields located near Gilmore City, Iowa. These plots were planted with rotations of corn (Zea mays L.) and soybean [Glycine max (L.) Merr.]. Corrugated plastic drainage pipe was installed at a depth of 1.06 m and a spacing of 7.6 m. The RZWQM2 model was calibrated manually by Qi et al. (2011). The parameters adjusted were mainly crop growth parameters, the soil root growth factor, and the hydraulic gradient. In the RZWQM2 model, crop growth parameters directly control the yield, biomass, and crop N uptake. Subsurface drainage is affected by crop growth because crop evapotranspiration is a function of leaf area index. Likewise NO 3 —N loss was tied to crop growth because N uptake is the major N sink in the model. The soil root growth factor determines crop water availability and uptake, which, in turn, indirectly influences drainage flow, N uptake, and N loss. Lateral water loss is estimated by the product of the hydraulic gradient and water table depth. Lateral flow loss and deep seepage are the two major factors that affect subsurface drainage flow. In our case, the deep seepage rate was set to zero. Parameters that were manually calibrated by Qi et al. (2011) and their ranges are listed in Table 1. In the present study, all parameters listed in Table 1 were recalibrated using the QPSO method against field observed data from CTRL1 treatment in Qi et al. (2011). The ranges for G2, G3, PHINT, and LFMAX represented the maximum and minimum values listed for all varieties in the DSSAT database. Soil root growth factors for each soil layer were computed using Eq. (11). In contrast to fixed values for wcg set by Ma et al. (2012), in this study, the range [1.0, 3.5] was used to calibrate it while running of QPSO algorithm until an appropriate solution was found, while z max was fixed at 100 cm. The initial values were randomly generated within the range of each parameter by the QPSO algorithm.

( SRGF ¼

1  wcg z 1  z max

z  15 cm z > 15 cm

ð11Þ

2.4. Objective function definition

Fig. 1. The procedure of QPSO algorithm.

Calibration of a model is always to be translated into the optimization of an objective function which contains the information of the calibrated model with respect to a certain problem. This objective function can then serve in programming for optimizing algorithms. The design of the objective function is therefore vital to an optimization problem, and a well suited objective function can lead to a well calibrated model. Given the wide range of field data values, the root mean square error (RMSE) with weighed coefficients was applied to design the objective function. Two group tests were carried out using the QPSO algorithm, one with three single object functions for each of three output components and another object function incorporating all the

75

M. Xi et al. / Computers and Electronics in Agriculture 113 (2015) 72–80 Table 1 Initialization range of crop parameters for QPSO algorithm. Crop

Parameter description

Initialization range

Corn

G2 G3 PHINT SRGF_wcg

Maximum possible number of kernels per plant Kernel filling rate during the linear grain filling stage and under optimum conditions (mg d1) Phylochron interval between successive leaf tip appearances Soil root growth factor parameter

[248, 990] [5, 16.5] [30, 50] [1.0, 3.5]

Soybean

LFMAX SRGF_wcg

Maximum leaf photosynthesis rate at 30 C, 350 vpm CO2, and high light (mg CO2 m2 s1) Soil root growth factor parameter

[0.8, 1.5] [1.0, 3.5]

Lateral flow

LHG

Lateral hydraulic gradient

[0, 104]

three output components. The ability of QPSO to calibrate RZWQM2 was first assessed with three single-objective functions representing grain yield (y), drainage flow (d), and NO 3 —N loss (nit), respectively. The three single-objective functions are similar and given as follows.

Objval y ¼

rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi i2 1 Xn h  obs wi Y 1;i  Y sim 1;i ðf ðÞÞ i¼1 n

ð12Þ

Objval d ¼

rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi i2 1 Xn h  obs wi Y 2;i  Y sim 2;i ðf ðÞÞ i¼1 n

ð13Þ

Objval nit ¼

rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi i2 1 Xn h  obs wi Y 3;i  Y sim 3;i ðf ðÞÞ i¼1 n

ð14Þ

where f() is a set of calibration parameters for RZWQM2, i denotes the year, ranging from 1 to 5, for the years 2005–2009, wi is the ratio of the parameter maximum over the five years to that in the particular year of interest:

wi ¼

 max Y obs i

ð15Þ

Y obs i

Objval_d is the optimization problem for drainage, Objval_nit is the optimization problem for NO 3 —N loss, Objval_y is the optimization problem for yield, obs obs Y obs 1;i ; Y 2;i ; Y 3;i are observed values, and  sim sim Y sim 1;i ; Y 2;i ; Y 3;i are values of yield, drainage flow and NO3 —N loss drawn from the output files of RZWQM2.

The calculated values for wi are given in Table 2. Secondly, in view of the complex and empirical nature of agriculture, users do not always be required to calibrate a model with one component, though it may be easier. Then one multipleobjective function which includes all the three output components was designed to calibrate the RZWQM2 model. This was accomplished by using the following equation.

rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi m i2 X 1 Xn h  obs sim Objval mul ¼ cj w j;i Y j;i  Y j;i ðf ðÞÞ i¼1 n j¼1

ð16Þ

Table 2 Weighted coefficients for each year of three single-objective functions. Year

Coefficient

Yield

Drainage

NO 3 —N

2005 2006 2007 2008 2009

w1 w2 w3 w4 w5

1.0 3.1 1.4 3.4 1.3

1.9 4.0 1.0 1.0 2.6

1.7 3.5 1.0 1.2 2.3

obs where n; i; wj;i ; Y sim j;i ; Y j;i ; f ðÞ are as defined for Eq. (15) and j is the jth component of the model, e.g., yield, and m has a value of 3 to denote the three components of yield, drainage flow and NO 3 —N loss.

To address the large difference in numerical magnitude between the observed values of yield, drainage and NO 3 —N loss (e.g., for 2005 yield was 10,200 kg ha1, but drainage water depth 25.8 cm) a weighed coefficient ck was introduced to balance the influences of the different components, so as to achieve equivalent parameter effects on the objective function. The mean observed values proportion method was used to calculate ck:

P5

ck ¼ Pi¼1 5

Y obs 1i

obs i¼1 Y ki

k ¼ 1; 2; 3

ð17Þ

Values of 1.0, 204.0, 153.4 were obtained for c1, c2, c3, respectively. For the multiple-objective function, if only c1 is kept (i.e., c1 = 1.0, c2 = 0, c3 = 0) in Eq. (16), then the multiple-objective function optimization translates into a single-objective function optimization (Eq. (12)). In the present study, twenty particles were used in the QPSO algorithm and the swarm was iterated 50 times. Coefficient b in the algorithm decreased from 1.0 to 0.5, which proved to be a better choice than most algorithm test functions (Sun et al., 2012) during the course of optimization. The best fitness values of swarm for each iteration were recorded. The parameters calibration procedure of RZWQM2 with the QPSO algorithm is shown in Fig. 2. The execution began with a population of parameters of RZWQM2 that were generated at random, according to the ranges given in Table 1. Each possible parameters solution was then written into the RZWQM2 input files and the performance of the model were tested. A new parameters population was then created from the existing population with QPSO algorithm and the process continued until the maximum number of iterations was achieved. 2.5. Model performance evaluation The performance of the QPSO-calibrated RZWQM2 model was compared with that of the manually-calibrated RZWQM2 model, using the same statistical analyses employed by Qi et al. (2011): percent bias (PBIAS), Nash–Sutcliffe efficiency (NSE), and the ratio of root mean squared error to standard error (RSR):

Pn  obs  Y sim  ð1 0 0Þ i i¼1 Y i PBIAS ¼ Pn  obs i¼1 Y i

ð18Þ

2 Pn  obs  Y sim i i¼1 Y i NSE ¼ 1  2 Pn  obs Y i¼1 Y i

ð19Þ

76

M. Xi et al. / Computers and Electronics in Agriculture 113 (2015) 72–80

Table 3 Calibrated parameters of RZWQM2 obtained through the QPSO algorithm. Crop

Parameter

Single-objective optimization

Multiple-objective optimization

Yield

Drainage

NO 3 —N

Corn

G2 G3 PHINT SRGF_wcg

319 12.58 30.10 3.45

944 6.68 48.71 2.86

602 8.15 49.99 1.00

284 14.11 48.66 1.41

Soybean

LFMAX SRGF_wcg

0.94 1.00

1.50 3.50

1.47 1.02

1.09 3.27

Lateral flow

LHG

2.07 105

3.85 107

4.25 1017

4.24 107

loss

Table 4 Yearly and 5-year (2005–2009) mean observed and simulated (RZWQM2 calibration with QPSO algorithm) annual yield, drainage flow and NO 3 —N loss for single-objective optimization problems.

rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2ffi Pn  obs sim  Yi i¼1 Y i RSR ¼ rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi  2ffi Pn obs Y  Y i i¼1

ð20Þ

where n is the number of data pairs, Y is the mean of the observed values, is the ith observed value, and Y obs i is the ith predicted value. Y sim i when NSE > 0.50, RSR < 0.70, and PBIAS is within ±25%, the performance of the model is considered to be ‘‘satisfactory’’ (Moriasi et al., 2007). The calibrated parameters were also applied to another treatment (CTRL2) in Qi et al. (2011) to validate the performance of QPSO. The validation CTRL2 data set is different from the calibration CTRL1 data set in terms of crops planted. In CTRL1 corn was planted odd years and soybean in even years, but in CTRL2 plots it was reversed with soybean in odd years and corn in even years.

3. Results and discussion The calibrated parameters of RZWQM2 for three singleobjective functions and one multiple-objective function with QPSO algorithm are presented in Table 3. For soil depths of 0.05, 0.15, 0.30, 0.45, 0.60, and 0.90 m, the soil root growth factors (SRGF) were calculated (Eq. (11)) for the single-objective optimization problems for yield [1, 1, 0.29, 0.13, 0.04, and 0, respectively], drain flow [1, 1, 0.36, 0.18, 0.07, and 0] and NO 3 —N losses [1, 1, 0.7, 0.55, 0.4, and 0.1] for corn, and [1, 1, 0.70, 0.55, 0.40, and 0.10], [1, 1, 0.29, 0.12, 0.04, and 0], [1, 1, 0.69, 0.54, 0.39, and 0.10] for soybean. For the same soil depths, the SRGF for the multiple-objective optimization problem were calculated to be [1, 1, 0.59, 0.43, 0.27, and 0.04] for corn and [1, 1, 0.31, 0.14, 0.05, and 0] for soybean.

The corresponding optimization results for yield, drainage and NO 3 —N loss from 2005 to 2009 are listed in Table 4. Over the 5-year period, simulated drainage depth and NO 3 —N losses were 28.0 cm and 40.3 kg N ha1, compared to measured values of 31.1 cm and 41.3 kg N ha1. Simulated and observed mean yield for corn was 8394 kg ha1 versus 8467 kg ha1 for corn, and were 3150 kg ha1 versus 3154 kg ha1 for soybean. Simulated annual corn and soybean yields (PBIAS = 0.7%, NSE = 0.98, and RSR = 0.20), drainage (PBIAS = 9.9%, NSE = 0.90, and RSR = 0.44), and NO 3 —N losses (PBIAS = 2.5%, NSE = 0.98, and RSR = 0.20) showed a close agreement with measured values. Thus all three evaluating indicators showed the model calibrated with the QPSO algorithm to perform in a ‘‘satisfactory’’ manner for all three problems. The evolutionary progress and variation of a single objective function during the progress of the QPSO algorithm through 50 iterations is shown for yield, drainage flow and NO 3 —N loss problems (Fig. 3A–C, respectively). All three function values showed fast convergence speed, such that after 20 iterations, the convergence rate slowed down as algorithm entered the local optimization stage. For the yield, the objective function value (Objval_y in Eq. (12)) decreased from 1380 to 595 after 30 iterations, and thereafter remained largely constant. Likewise, the objective function value for drainage (Objval_d in Eq. (13)) dropped from 7.24 to 6.47 after 20 iterations for the drainage flow simulation, and thereafter remained largely constant. Similarly, the objective function value for NO 3 —N loss (Objval_n in Eq. (14)) showed a similar trend slowing after 20 iterations, and remained constant thereafter. Progress in the values of PBIAS, NSE and RSR indicators during the optimization procedure are shown in Figs. 3D–F, respectively. In general, the performance of the model for simulating yield and NO 3 —N loss was improved by the QPSO algorithm based on single-objective functions. For example, the magnitude of PBIAS dropped from 2% to 0.6% for yield and from 7.5% to 3.5% for NO 3 —N loss. However, the statistical indicators do not suggest an obviously improvement of simulation in subsurface drainage flow,

M. Xi et al. / Computers and Electronics in Agriculture 113 (2015) 72–80

Fig. 2. Calibration procedure for RZWQM2 with QPSO algorithm.

though the objective function value of drainage showed an evident decrease during the iteration. This may be attributed to the fact that the search for the index of the LHG was in the range of [1, 4], and in the first iteration, it is easy for one of 20 particles in the swarm to fall into the range within which the optimal solution is located. However, if the search for LHG was in a range of [0, 0.0001], the particles would need more iterations (about 200) to map the optimal solution, thus significantly increase the necessary running time. Therefore, we suggest searching for the index for LGH when optimizing for the drainage problem. Although the QPSO method did not improve the simulation results in comparison to the manual calibration method in this study, QPSO is still useful because it is a subjective approach and saves man power. 3.1. Multiple-objective optimizations While the QPSO algorithm displayed the ability to satisfactorily calibrate RZWQM2 in single-objective optimizations, users most often need to check most or all system outputs, not just one with the corresponding field measurements. When there are many outputs, the use of a multiple-objective optimization may be

77

useful. Therefore, the QPSO algorithm was applied to a multipleobjective optimization problem for RZWQM2 parameter calibration, a more complex process than for a single-objective optimization. In this group test, the QPSO algorithm had the same parameters settings as for the single-objective optimization, i.e., 20 particles, 50 iterations and b with a value decreasing from 1.0 to 0.5. Optimization results for yield, drainage flow and NO 3 —N loss (Table 5) show simulations to be generally satisfactory for loss. (11.2% < PBIAS < 1.5%, yield, drainage, and NO 3 —N NSE > 0.80, and RSR 6 0.41 (Table 5). The statistical indicators were comparable with the results by manual calibration in Qi et al. (2011). Corresponding calibration parameters of RZWQM2 are listed in ‘‘multiple objective’’ column of Table 3. Over the 5-year period, simulated mean yield was 7899 kg ha1 and 3172 kg ha1, respectively, for corn and soybean, comparing to the observed value of 8470 kg ha1 and 3150 kg ha1. Simulated drainage depth and NO losses were 27.6 cm and 3 –N 41.7 kg N ha1, compared to measured values of 31.1 cm, and 41.3 kg N ha1. Values of PBIAS, NSE, and RSR for yield were 5.2%, 0.85, and 0.38, respectively, while those for NO 3 —N loss were 1.5%, 0.98, and 0.13 (Table 5), in both cases close to those obtained using the manually calibrated model (Qi et al., 2011). However, the QPSO-assisted model’s performance for drainage was not as good as that of the manually-calibrated model. This could be attributed to differences in soil root growth factors (SRGF): under manual calibration the SRGF of each soil layer was adjusted by trial and error, whereas under the QPSO method, the SRGF were computed using Eq. (11) by calibrating wcg within the range of [1.0, 3.5]. The QPSO algorithm resulted in higher values of SRGF, which led to higher water uptake. For example, the QPSO algorithm calibrated SRGF values were 1, 1, 0.59, 0.43, 0.27, and 0.04 for corn while the values obtained by manual calibration were 1, 1, 0.3, 0.07, 0.07, and 0.01. However, this did not result in higher corn crop water transpiration. Because the PHINT parameter optimized with the QPSO algorithm was higher than that obtained under manual calibration, both the number of leaves and leaf area index were 12.5% and 13.5% less, respectively than that for manual calibration for corn. These low leaf area indices resulted in a 6.7% lower plant water uptake, which explained the lower simulated drainage for the manual calibration. Nonetheless, the SRGF values obtained by the QPSO algorithm were more reasonable than those smaller values obtained by manual calibration (Ma et al., 2012). Fig. 4A illustrated the evolutionary progress of multipleobjective function during optimization with the QPSO algorithm. The objective value continually declined from 6091 to 3134. Sometime it would stop at a local optimized position, but jumped out quickly and then kept on searching for the global best solution. For example, during iteration steps of 7 and 11, the swarm was trapped in a local optimal area, but at iteration step 12, the swarm escaped from the area and move to a better solution. This repeatedly happened during the course of optimization. In Fig. 4B–D, the value of PBIAS, NSE, and RSR varied over the course of evolution, but always were finding towards a more reasonable value. For a certain indicator, we did not get the best value at the 50th iteration. For example, the PBIAS of yield was better with value 1.03% at 30th iteration while it is 5.23% at 50th. This was due to the optimization problem which considered all three components including yield, drainage flow and NO3–N loss as one multiple-objective function. During the course of evolution, the three aspects will be considered and when the PBIAS value of yield turned to be worse at some iterations, the PBIAS value of the other components will be improved, and the whole model tended to be better and better. We could find that the PBIAS of water drainage and NO3–N loss were 9.3% and 7.6%, respectively at the 30th iteration, and these two values are 11.22% and 1.50% at 50th iteration respectively.

M. Xi et al. / Computers and Electronics in Agriculture 113 (2015) 72–80

1200 1000 800 600

A

7.2 7.0 6.8 6.6

B

10

20

30

40

0

50

Objective fuction value of drainage

0.05 0.00 -0.05 -0.10

D

-0.15

10

20

30

40

20

30

40

0.9 0.8 0.7 0.6

50

10 8 6 4

C

E 0

0

50

1.0

0.5 0

10

Objective fuction value of NO3 -N loss

0

12

2

6.4

400

Objective fuction value of yield

7.4

Objective fuction value of NO3 -N loss

1400

Objective fuction value of drainage

Objective fuction value of yield

78

10

20

30

40

50

0.8 Yield Drainage NO 3-N loss

0.6 0.4 0.2

F 0.0

10

Iterations

20

30

40

50

0

10

Iterations

20

30

40

50

Iterations

1 Fig. 3. Evolutionary progress of (A) yield (kg ha1), (B) drainage water (cm yr1), (C) NO ), and for each output parameter evolution of (D) PBIAS, (E) NSE, and 3 —N loss (kg ha (F) RSR values for three single optimization problems addressed with a QPSO algorithm.

0.8

1.0

7000 0.1

6000

0.9 0.6

3000 2000

0.0

-0.1

0.8

RSR value

4000

NSE value

5000

PBIAS value

Objective function values

Table 5 Mean observed and manually (Qi et al., 2011) or QPSO-assisted simulated crop yield, drainage depth and NO 3 —N loss in CTRL1 treatment approached as a multiple-objective optimization problem.

0.7 0.6 Yield Drainage NO3-N loss

0.5 1000

A 0

10

B

-0.2 20

30

Iterations

40

50

0

10

C

0.4 20

30

Iterations

40

50

0

10

20

30

Iterations

40

50

0.4

0.2

D

0.0 0

10

20

30

40

50

Iterations

Fig. 4. Evolutionary progress of (A) objective function values, (B) PBIAS, (C) NSE and (D) RSR in a multiple-objective optimization problem with QPSO algorithm.

M. Xi et al. / Computers and Electronics in Agriculture 113 (2015) 72–80

79

Table 6 Validation of the performance of QPSO-assisted calibrated parameters using data set from the CTRL2 treatment.

For the objective function value the same conclusion could be reached since the value was 4087 at 30th iteration but only 3134 at 50th iteration. So in the test, QPSO algorithm was always searching for a good set of parameters to balance the three aspects. Compared with single-objective optimization problems, multiple-objective problems converged more slowly. After 30 iterations, all the single-objective problems entered into a global optima area, while for the multiple-objective problem 46 iterations were required. This suggests that more iterations are needed for the QPSO algorithm when a more complex objective function is established. One should note that the QPSO algorithm calibrated G2 and G3 values were very different from that obtained by manual calibration (Qi et al., 2011; Ma et al., 2012). The maximum kernel number (G2) found by QPSO was 284, while in most manual calibration cases it ranges from 700 to 1000. The low G2 was compensated by a higher kernel filling rate (G3) in order to obtain a good simulation of corn yield. The low G2 obtained by QPSO algorithm can be attributed to the wide range that was set for this parameter. The range of G2 was selected using the minimum and maximum values found in the DSSAT database for all corn varieties. Therefore, experience is needed to set a reasonable range for each parameter for a given problem. Figs. 2 and 4 both show that the calibrated parameters obtained by the first iteration resulted in a satisfactory or near satisfactory model performance. This can be explained by the fact that the probability of one randomized particle in the swarm falling into an area of satisfactory solutions is high. It more obviously happened for the single-objective problem than the multipleobjective function. The total time for conducting 50 iteration using 20 particles for the 5 years of data was about 50 h on an IBM X230 laptop computer.

4. Validation The validation results showed that the performance of QPSO calibrated parameters was comparable to the manual calibration (Table 6). The yield was overestimated more by the QPSO calibrated parameters than those manual calibrated by Qi et al. (2011), but the simulations in drainage and NO 3 –N loss were slightly better. The overestimation in corn and soybean yield suggested that QPSO calibrated G2, G3, and LFMAX was slightly high. However, the simulation in NO 3 –N loss using QPSO calibrate parameter set was much better than the manual calibration, in particular in 2005 and 2006. The differences in predicted drainage between QPSO and manual calibrated parameters were minor.

5. Conclusions The present study explored the application of a quantum-behaved particle swarm optimization (QPSO) algorithm in the automatic calibration of the RZWQM2 model. RZWQM2 was automatically calibrated for a field sub-dataset collected from a subsurface drained corn-soybean rotation in Iowa, which included yield, drainage flow and NO3–N loss with two different functions. One was focused on the three components separately (singleobjective) and the other considered all three aspects (multipleobjective). The calibrated RZWQM2 models were satisfactory with percent bias (PBIAS) values within ±25%, Nash–Sutcliffe efficiency (NSE) > 0.50, and a ratio of root mean square error to standard error (RSR) < 0.70. The results demonstrate that the QPSO algorithm can serve as a valuable tool in automatically calibrating the RZWQM2 model. More iterations are needed when a more complicated objective function is being considered, which would require longer processing time. User experience is necessary for setting appropriate ranges for the model parameters. Possible future studies include (i) the design of a paralleled QPSO algorithm to reduce the calibration time of RZWQM2 according to the parallel characteristics of the QPSO algorithm; and (ii) testing the use of the QPSO algorithm to calibrate other agriculture models. Acknowledgements The research work was supported by the Jiangsu Overseas Research & Training Program for University Prominent Young & Middle-aged Teachers and Presidents, and partly supported by the Natural Science Foundation of Jiangsu Province, China (Project Number: BK2010143), by the National Natural Science Foundation of China (Project Numbers: 61170119 and 61105128), and by the Program for New Century Excellent Talents in University (Project Number: NCET-11-0660), as well as the Qing Lan Project of Jiangsu. References Ahuja, L.R., Rojas, K.W., Hanson, J.D., Shaffer, M.J., Ma, L., 2000. Root Zone Water Quality Model: Modeling Management Effects on Water Quality and Crop Production. Water Resources Publications, Highlands Ranch, Colo. Angeline, P.J., 1998. Evolutionary optimization versus particle swarm optimization: philosophy and performance differences. p. 601–610. In: Porto, V.W., Saravanan, N., Waagen, D., Eiben, A.E. (Eds.), Evolutionary Programming VII. Proceedings of the Seventh International Conference, EP98, San Diego, California, USA, March 25–27, Springer, Berlin, Heidelberg. DOI:http://dx.doi. org/10.1007/BFb0040811. Calmon, M.A., Jones, J.W., Shinde, D., Specht, J.E., 1999. Estimating parameters for soil water balance model using adaptive simulated annealing. Appl. Eng. Agric. 15 (6), 703–713. http://dx.doi.org/10.13031/2013.5841.

80

M. Xi et al. / Computers and Electronics in Agriculture 113 (2015) 72–80

Clerc, M., Kennedy, J., 2002. The particle swarm-explosion, stability and convergence in a multidimensional complex space. IEEE Trans. Evol. Comput. 6 (2), 58–73. http://dx.doi.org/10.1109/4235.985692. Coelho, L.S., Mariani, V.C., 2008. Particle swarm approach based on quantum mechanics and harmonic oscillator potential well for economic load dispatch with valve-point effects. Energy Convers. Manage. 49 (11), 3080–3085. http:// dx.doi.org/10.1016/j.enconman.2008.06.009. Fang, Q.X., Green, T.R., Ma, L., Erskine, R.H., 2010. Optimizing soil hydraulic parameters in RZWQM2 under fallow conditions. Soil Sci. Soc. Am. J. 74 (6), 1897–1913. http://dx.doi.org/10.2136/sssaj2009.0380. Gao, F., Tong, H.-Q., 2006. Parameter estimation for chaotic system based on particle swarm optimization. Acta Phys. Sin. 55 (2), 577–582. Jones, J.W., Hoogenboom, G., Porter, C.H., Boote, K.J., Batchelor, W.D., Hunt, L.A., Wilkens, P.W., Singh, U., Gijisman, A.J., Ritchie, J.T., 2003. The DSSAT cropping system model. Eur. J. Agron. 18 (3–4), 235–265. http://dx.doi.org/10.1016/ S1161-0301(02)00107-7. Kennedy, J., Eberhart, R.C., 1995. Particle swarm optimization. In: Proceedings of the IEEE International Conference on Neural Networks, vol. 4, Nov. 27–Dec. 1, 1995, Piscataway, NJ, pp. 1942–1948. DOI:http://dx.doi.org/10.1109/ICNN.1995. 488968. Kuok, K.K., Chan, C.P., 2012. Particle swarm optimization for calibrating and optimizing Xinanjiang model parameters. Int. J. Adv. Comput. Sci. Appl. 3 (9), 115–123. http://dx.doi.org/10.14569/IJACSA.2012.030917. Lei, X., Fu, A., 2008. Two-dimensional maximum entropy image segmentation method based on quantum-behaved particle swarm optimization algorithm. In: Fourth International Conference on Natural Computation, vol. 3, 18–20 October 2008, Jinan, PRC, pp. 692–696. DOI:http://dx.doi.org/10.1109/ICNC.2008.822. Ma, L., Hoogenboom, G., Ahuja, L.R., Nielsen, D.C., Ascough II, J.C., 2005. Development and evaluation of the RZWQM–CROPGRO hybrid model for soybean production. Agron. J. 97 (4), 1172–1182. http://dx.doi.org/10.2134/ agronj2003.0314. Ma, L., Hoogenboom, G., Ahuja, L.R., Ascough II, J.C., Saseendran, S.A., 2006. Evaluation of the RZWQM–CERES-Maize hybrid model for maize production. Agric. Syst. 87 (3), 274–295. http://dx.doi.org/10.1016/j.agsy.2005.02.001. Ma, L., Trout, T.J., Ahuja, L.R., Bausch, W.C., Saseendran, S.A., Malone, R.W., Nielsen, D.C., 2012. Calibrating RZWQM2 model for maize response to deficit irrigation. Agric. Water Manag. 103, 140–149. http://dx.doi.org/10.1016/ j.agwat.2011.11.005. Malone, R.W., Jaynes, D.B., Ma, L., Nolan, B.T., Meek, D.W., Karlen, D.L., 2010. Soiltest N recommendations augmented with PEST-optimized RZWQM simulations. J. Environ. Qual. 39 (5), 1711–1723. http://dx.doi.org/ 10.2134/jeq2009.0425. Mikki, S., Kishk, A.A., 2006. Infinitesimal dipole model for dielectric resonator antennas using the QPSO algorithm. In: Antennas and Propagation Society International Symposium 2006, IEEE, 9–14 July 2006, Albuquerque, NM, pp. 3285–3288. DOI:http://dx.doi.org/10.1109/APS.2006.1711314. Moriasi, D.N., Arnold, J.G., Van Liew, M.W., Binger, R.L., Harmel, R.D., Veith, T.L., 2007. Model evaluation guidelines for systematic qualification of accuracy in watershed simulations. Trans. ASABE 50 (3), 885–900. http://dx.doi.org/ 10.13031/2013.23153. Nolan, B.T., Puckett, L.J., Ma, L., Green, C.T., Bayless, E.R., Malone, R.W., 2010. Predicting unsaturated zone nitrogen mass balances in agricultural settings of the United States. J. Environ. Qual. 39 (3), 1051–1065. http://dx.doi.org/ 10.2134/jeq2009.0310. Pabico, J.P., Hoogenboom, G., McClendon, R.W., 1999. Determination of cultivar coefficients of crop models using a genetic algorithm: a conceptual framework. Trans. ASAE 42 (1), 223–232. http://dx.doi.org/10.13031/2013.13199.

Poli, R., Kennedy, J., Blackwell, T., 2007. Particle swarm optimization an overview. Swarm Intell. 1 (1), 33–57. http://dx.doi.org/10.1007/s11721-007-0002-0. Qi, Z., Helmers, M.J., Malone, R.W., Thorp, K.R., 2011. Simulating long-term impacts of winter rye cover crop on hydrological cycling and nitrogen dynamics for a corn–soybean crop system. Trans. ASABE 54 (5), 1575–1588. http://dx.doi.org/ 10.13031/2013.39836. Sabat, S.L., Coelho, L.S., Abraham, A., 2009. MESFET DC model parameter extraction using quantum particle swarm optimization. Microelectron. Reliab. 49 (6), 660– 666. http://dx.doi.org/10.1016/j.microrel.2009.03.005. Shaffer, M.J., Rojas, K.W., Decoursey, D.G., Hebson, C.S., 2000. Chapter 5: nutrient chemistry processes – OMNI. In: Ahuja, L.R., Rojas, K.W., Hanson, J.D., Shaffer, M.J., Ma, L. (Eds.), Root Zone Water Quality Model: Modeling Management Effects on Water Quality and Crop Production. Water Resources Publications, Highlands Ranch, CO, pp. 119–144. Shi, Y., Eberhart, R.C., 1998. Parameter selection in particle swarm optimization. In: Evolutionary Programming VII. Proceedings of the Seventh International Conference, EP98 (San Diego, CA, March 25–27, 1998), Springer, Berlin, Heidelberg, pp. 591–600. DOI:http://dx.doi.org/10.1007/BFb0040810. Shuttleworth, W.J., Wallace, J.S., 1985. Evaporation from sparse crops – an energy combination theory. Q. J. R. Meteorol. Soc. 111 (469), 839–855. http:// dx.doi.org/10.1002/qj.49711146910. Sun, J., Feng, B., Xu, W., 2004a. Particle swarm optimization with particles having quantum behavior. In: Proceedings of Congress on Evolutionary Computation, 2004, CEC2004 (Portland OR, 20–23 June 2004), vol. 1, pp. 326–331. DOI:http:// dx.doi.org/10.1109/CEC.2004.1330875. Sun, J., Feng, B., Xu, W., 2004b. A global search strategy of quantum-behaved particle swarm optimization. In: Proceedings of 2004 IEEE Conference on Cybernetics and Intelligent Systems (Singapore, 1–3 Dec. 2004), vol. 1, pp. 111– 116. DOI:http://dx.doi.org/10.1109/ICCIS.2004.1460396. Sun, J., Liu, J., Xu, W., 2006. QPSO-based QoS multicast routing algorithm. In: Wang, T.-D., Li, X., Chen, S.-H., Wang, X., Abbass, H., Iba, H., Chen, G.-L., Yao, X. (Eds.), Simulated Evolution and Learning. Proceedings 6th International Conference, SEAL 2006 (Hefei, China, October 15–18, 2006), Springer, Berlin, Heidelberg, pp. 261–268. DOI:http://dx.doi.org/10.1007/11903697_34. Sun, J., Fang, W., Wu, X., Vasile, P., Xu, W., 2012. Quantum-behaved particle swarm optimization analysis of individual particle behavior and parameter selection. Evol. Comput. 20 (3), 349–393. http://dx.doi.org/10.1162/EVCO_a_00049. Tang, L., Xue, F., 2008. Using data to design fuzzy system based on quantumbehaved particle swarm optimization. In: International Conference on Machine Learning and Cybernetics (Kunming, PRC, 12–15 July 2008), vol. 1, pp. 624–628. DOI:http://dx.doi.org/10.1109/ICMLC.2008.4620479. Thorp, K.R., Malone, R.W., Jaynes, D.B., 2007. Simulating long-term effects of nitrogen fertilizer application rates on corn yield and nitrogen dynamics. Trans. ASABE 50 (4), 1287–1303. http://dx.doi.org/10.13031/2013.23640. Xi, M., Sun, J., Xu, W., 2006. Quantum-behaved particle swarm optimization for design H1 structure specified controllers. In: Wenbo, X., Guangming, W. (Eds.), International Symposium on Distributed Computing and Applications to Business, Engineering and Science (DCABES), vol. 2, pp. 1016–1019 (seen 2 March 2014 at ). Xi, M., Sun, J., Xu, W., 2007. Parameter optimization of PID controller based on quantum-behaved particle swarm optimization. In: Complex Systems and Applications – Modeling, Control and Simulations. DCDIS Series B, vol. 14, pp. 603–607 (seen 2 March 2014 at .