Forecasting effluent quality of an industry wastewater treatment plant by evolutionary grey dynamic model

(1)

Contents lists available atScienceDirect

Resources, Conservation and Recycling

j o u r n a l h o m e p a g e :w w w . e l s e v i e r . c o m / l o c a t e / r e s c o n r e c

Forecasting efﬂuent quality of an industry wastewater treatment plant by

evolutionary grey dynamic model

Ho-Wen Chen

a,∗

_{, Ruey-Fang Yu}

b

_{, Shu-Kuang Ning}

c

_{, Hsin-Chih Huang}

a

a_{Department of Environmental Engineering and Management, Chao-Yang University of Technology, Taichung 413, Taiwan, ROC} b_{Department of Safety, Health and Environmental Engineering, National United University, Miao-Li 360, Taiwan, ROC} c_{Department of Civil and Environmental Engineering, National University of Kaohsiung, Kaohsiung, Taiwan, ROC}

a r t i c l e i n f o

Article history: Received 6 January 2009

Received in revised form 6 August 2009 Accepted 12 August 2009

Keywords: Grey systems theory Genetic algorithms Wastewater treatment plant Forecasting

Monte Carlo simulation

a b s t r a c t

The application of a prediction model is a commendable exercise to evaluate a facility’s performance and achieve better quality control in the operation of wastewater treatment plants. This paper proposes a model which integrates grey dynamic modeling and genetic algorithm to predict accurately the efflu-ent quality of an industrial wastewater treatmefflu-ent plant located in southern Taiwan. Model parameters, variables and structures are determined endogenously to minimize errors between observed and pre-dicted values. Modeling feasibility has been proved by using data compared with Monte Carlo simulation and artificial neural network approaches. The results show that the prediction of our proposed model is sensitive to the joint effect of suspended solids and F/M.

1. Introduction

In order to meet increasingly stricter government regulation and standard operating procedure, a wastewater treatment plant (WTP) operator is obliged to collect operation data for use in facilities per-formance evaluation and for better control of an operation for stable efﬂuent quality. To maintain a stable efﬂuent quality, WTP oper-ators are required to adhere to standard operating procedures for routine control and emergency response. Once a well-trained oper-ator leaves the position or otherwise becomes absent, replacement operators often face tough challenges in performing procedures in accordance to the stricter standards requirements. In response to these concerns, many researchers have engaged in developing use-ful prediction models to facilitate operators’ tasks (Lee et al., 2005; Gernaey et al., 2004). The efforts based on white-box, black-box, stochastic grey-box and hybrid techniques have made many contri-butions to conquering the complex nature of wastewater treatment (Belanchea et al., 2000; Hanbay et al., 2008; Mjalli et al., 2007; Tenno and Uronen, 1995). By and large, reliability is the major con-cern for operators handling WTP operations. It is axiomatic that the more accurate the model is, the larger the data collection efforts must be, assuming that all other things are equal. So the opera-tion of a large WTP entails massive efforts in computaopera-tion and data

∗ Corresponding author. Tel.: +886 4 23323000x4395; fax: +886 4 23742365. E-mail address:[email protected](H.-W. Chen).

collection for building complex prediction models, such as neural network and time series models (Van Dongen and Geuens, 1998; Berthouex and Box, 1996). Given the fact that data collection is time-consuming and expensive, an operator will invariably exer-cise discretion as to which quality parameters to utilize and how much data to collect.

For some developing countries, like Indonesia, Philippine, India, China, there are many small- and mid-scale municipal wastewater treatment plants (MWWTP) built to resolve serious problems in water quality protection and water resource conservation. Accord-ing to statistics byWei (2000), the number of MWWTPs in China will have rapidly increased by over 1000 by the year 2010, when most of WTPs in China will be small- to mid-scale. Retaining a sta-ble efﬂuent quality for these small- and mid-scale WTPs might be a strict challenge for inexperienced operators, especially when mon-itoring data is incomplete. In fact, most medium- and small-size WTPs cannot afford to expend massive efforts to pursue a perfect model. Therefore a host of simpliﬁed methods has been proposed for use in situations lacking data and where doubtful model esti-mation parameters exist (Sophie and Adeline, 2007; Marsili-Libelli, 2004; Belanche et al., 1999).

For instance, since its use in conjunction with a difference equa-tion byDang to simulate energy systems in 1982, grey dynamic modeling (GM) has been applied to environments with poor, incomplete, or uncertain information (Hsu and Chen, 2003; Chen and Chang, 2000). One reason for these applications is that a GM model demands no explicit formulation for physiochemical

(2)

Fig. 1. The schematic of wastewater treatment process.

properties or socioeconomic behavior (Leephakpreeda, 2005; Lin and Lin, 2001; Mao and Chirwa, 2006; Zhou et al., 2006).Akay and Atak (2007) applied a grey model with a rolling mecha-nism approach to predict Turkey’s total and industrial electricity consumption.Wordena et al. (2007)present grey-box and black-box-based models to forecast the slipping time of a landslide and to successfully predict a landslide. Studies have also been done for developing hybrid soft forecasting techniques (Chang and Tsai, 2008). These include the studies ofOlson and Wu (2006)to inte-grate fuzzy and grey dynamic models (Wang, 2004) andHe and Hwang (2007)who use genetic algorithms to optimize the param-eters of GM modeling (Wang and Hsu, 2008; Tarng et al., 2002).

As with most model-based mathematical modeling, the expe-rience of model builders is essential to the determination of instrumental variables and parameters. This paper uses genetic algorithm (GA) to establish a data-based GM model for a given monitoring database. Variables and parameters modeled were determined after a series of evolutions was calculated by genetic algorithms. In addition, Monte Carlo simulation was used to quan-tify the sensitivity stemming from each decision variable. 2. Materials and methodology

This study’s researchers used an industrial wastewater treat-ment plant as a case study to build an evolutionary grey prediction model for forecasting the effluent concentration of chemical oxy-gen demand (COD). This WTP is an activated sludge process that handles the wastewater from a steel factory located in southern Taiwan. Its treatment train consists of equalization, aeration and solid–liquid separation, as indicated inFig. 1. The raw wastewater is composed of production line effluent, laboratory discharges and domestic sewage from within the factory. Because inflow quality varies over time, it is a real challenge to operate and control the plant correctly to meet the effluent requirements. Four water qual-ity items and operation variables are exhibited in a defined system (Fig. 1). Monitoring points for these variables are also shown in Fig. 1. With regard to the operating parameters and effluent quali-ties, 68 daily monitoring data were collected from January 2000 to June 2001, using the concentration of COD in the treatment effluent as the indicator for assessing the accuracy of a prediction model. 2.1. The principle of grey dynamic model

Grey systems theory posits that outcomes from a dynamic system will encompass some information that implicates the behaviors of the system. Therefore, instead of analyzing the char-acteristics of the grey systems directly, the grey systems theory

exploits the accumulated generating operation (AGO) technique to outline system behaviors, because external and latent behaviors will become more apparent due to the decrease of the random intensity after the AGO practice; that is, the AGO practice will reduce the white noise in a statistical sense.

If the original series of data with n samples is X(0) = (x(0)_(1),

x(0)_{(2), . . ., x}(0)_{(n)), then the output of this AGO is expressed as a}

ﬁrst-order accumulated generating sequence, X(1)_{= (x}(1)_{(1), x}(1)_(2),

. . ., x(1)_(n_{− 1), x}(1)_{(n)), where x}(1)_{(k) corresponds to the system}

out-put at time period k. The generation of higher order accumulated generating sequences may be obtained after an iterative procedure by following the same logic; that is, the r-order AGO series X(r)_can

be deﬁned as X(r)_(n)₌

n

k=1x(r−1)(k) for k = 1, 2, . . ., n.

GM(p, n) is the general form of a GM model, as shown as Eq.(1), in which p stands for the order of differential equations and n is the number of variables used for building the GM model.

p

i=0 ai dp−iX_i(1) dtp−i = n

i=2 biXi(1)+1 (1)

Eq.(2)is the speciﬁcation of GM(1, n), a stylized GM model, that is used to predict the outcome of a ﬁrst-order system with n variables, and Eqs.(2) and (3)are the discrete form of the model. a1 dx(1)₁ (k) dt + a0x (1) 1 (k)= n

i=2 bix_i(1)(k) (2) x(0)₁ (k)+ aZ(k) = n

i=2 bix(1)i (k) (3) Z(k)=1 2[x (1) 1 (k− 1) + x (1) 1 (k)] (4)

In this formulation, X(1)_{(k) is the accumulated value of X}(0)_from

time 1 to time k; parameters a and bjare obtained by using

least-squares method when the GM(1, n) model is viewed as a regression problem. The estimator, ˇ, for these parameters is listed in Eqs. (5)–(7). Z(k) is the background of the differential equation. ˇ= (a, b2, b3,· · ·, bn)= (BTB)−1BTyn (5) B=

⎡

⎢

⎣

−0.5(x(1) 1 (1)+ x (1) 1 (2)) x (1) 2 (2) · · · x (1) n (2) −0.5(x(1) 1 (2)+ x (1) 1 (3)) x (1) 2 (3) · · · x (1) n (3) . . . ... ... −0.5(x(1) 1 (k− 1) + x (1) 1 (k)) x (1) 2 (k) · · · x (1) n (k)

⎤

⎥

⎦

(6)

(3)

Fig. 2. The ﬂowchart of the establishment of evolutionary GM model.

yn= (x(0)(2), x(0)(3),· · ·x(0)(k)) T

(7) While ˇ is determined by Eqs.(5)–(7), difference equation of GM model can be resolved as Eq.(8)and(9). The solutions include two parts. The fitting output of ˆx(0)₁ (k+ 1) at time k + 1 is first obtained by Eq.(5), and then Eq.(8)is used as a manipulation basis to perform an inverse accumulated generating operation (i.e., denoted as [AGO]−1 or IAGO) by which the final result, ˆx(0)₁ (k+ 1) yields in Eq.(9)

ˆx(1)₁ (k+ 1) = (x0 1(1)− n

i=2 bix(1)_i (k) a )e −ak₊ n

i=2 bix(1)_i (k) a (8) ˆx(0)₁ (k_{+ 1) = ˆx}(0)₁ (k_{+ 1) − ˆx}(0)₁ (k) (9)

2.2. Evolutionary grey dynamic model

Fig. 2illustrates the analytical procedure for building an evolu-tionary GM(1, n) model, in which there are five main steps to be followed. In the first step, a first-order accumulated series is for-mulated by AGO practice. The goal of step 2 is to build an integer nonlinear programming instead of regression analysis used in tradi-tional GM models, as well as to decide the variables and parameters for modeling GM. Steps 3 and 4 are performed to resolve the inte-ger nonlinear programming by GA owing to its ability in searching for a global solution. In Step 5, IAGO is performed to generate the predicted values.

2.2.1. Objective function

To optimize the GM(1, n) prediction model by GA, root mean square error (RMSE), Eq.(10), is used as the objective function for a nonlinear integer programming with parameters and variable as decisions variables. min RMSE=

k (ˆx(0)₁ (k)− x(0) 1 (k)) 2 k− 1 (10)

2.2.2. Constraints for GM model

This paper uses binary variables in the optimization process to decide whether certain variables should be included in the GM model. Constraints are aided by Eqs. (11)–(13), in which giand

w_n−1+u+v are integer variables that decide whether the variables

in Eq.(11)should be added into a grey differential equation. Eq. (12)is to avoid over ﬁtting, in which Vmaxis the maximum number

of variables allowed in the GM model. dx(0)₁ dt + ax (1) 1 = n

i=2 gibixi(1)+1+ n

u=2 n

v=u+1 wn−1+u+vbn−1+u+vx(1)u xv(1) (11) n

i=2 gi+ n

u=2 n

v=u+1 wn−1+u+v≤ Vmax (12) gi, wn−1+u+v∈ {0, 1}

∀

i, n, u,

v

(13)

2.2.3. Accuracy improvement constraints

The background coefﬁcient, Z(k), expressed by Eq.(4), is an inﬂu-ential factor of the GM model. For better results, this paper applies a procedure with a three-point weighted average in the estimation. Z(1)_(k_{+ 1) = z}

1x(1)1 (k− 1) + z2x1(1)(k)+ z3x(1)1 (k+ 1) (14)

z1+ z2+ z3= 1 (15)

0≤ z1, z2, z3≤ 1 (16)

where z1, z2, and z3represent the corresponding weights of data in

different periods. 2.3. Solution procedure

GA implementation requires all decision variables and param-eters to be coded as chromosomes, a group of binary strings, and executes a series of evaluations and comparisons, including the three classic genetic operators in reproduction, crossover, and mutation, until the global optimal solution is found. Fitness func-tion, Eq.(9), is used as a criterion to investigate the improvement of the optimization process over generations. In this paper mutation rate is equal to 0.2 and crossover rate is equal to 0.35.

(4)

Table 1

Coefﬁcients and their statistical characteristics in evolutionary GM prediction model.

Variables Coefﬁcient Standard deviation P-value

COD a = 0.11 0.02 0.00000431* SSa _b 2=−348.81 73.08 0.00001235* HRTb _b 3= 68.82 15.54 0.00004164* R%c _b 4= 0e – – F/Md _b 5= 0e – – SS× F/M b6= 20.29 6.59 0.00316576* SS× R% b7= 36.41 7.29 0.00000558* SS× HRT b8= 0e – – HRT× R% b9=−7.73 1.51 0.00000362* HRT× F/M b10=−2.17 1.21 0.07833491 R%× F/M b11= 0e – – Note: *P-value≤ 0.01. a_{Suspended solids.} b_{Hydraulic retention time.} c _{Percentage of return sludge.} d _{Food-to-microorganism ratio.} e_{Variable that does not be selected by GA.}

2.4. Monte Carlo simulation method for sensitivity analysis

To assess the random errors and examine the stability of the proposed prediction model, Monte Carlo simulation was used with varying values for input variables. For ease of analysis and discus-sion, the monitoring items were grouped into operating variables and status variables according to their capacity. Based on mon-itoring data, probability distributions of, and cross correlation coefﬁcient among, these variables were estimated for use in gen-erating input data for simulation.

3. Results and discussion 3.1. Variable selection

After a series of evolutions by GA, six variables, shown inTable 2, were selected for model building. The coefﬁcients and their sta-tistical characteristics were presented inTable 1. Except for the

interactive terms, HRT and F/M, these variables are statistically signiﬁcant. R% and F/M did not exert their effects directly, but indi-rectly by joining with other variables on efﬂuent quality prediction. The variations of observed, predicted values and the mean absolute percentage error (MAPE) inFig. 3show that the evolutionary GM model achieves an accurate prediction with 19% error on average. The largest error happened on the 15th day, corresponding to an extremely high observation value.

For comparison purposes, stepwise regression analysis was applied to the data set.Table 2shows that all variables determined by stepwise regression are statistically significant if the number of predictor variables is six or less.Fig. 4presents the result of stepwise-regression-based model, in which HRT, R%, SS× HRT, and SS× R% were selected for modeling. Based on the comparison of Fig. 3withFig. 4, it appears that step-wise regression performs bet-ter than GM model, specifically at peaks. However, the evolutionary GM model has presented its ability in a small-value prediction. Clearly, the findings indicate that a step-wise, regression-based

Table 2

Selected variables for building GM model by stepwise regression model and genetic algorithm.

Numbers of variable

Average absolute error (%)

Adjusted-R2 _{Selected variables and their signiﬁcances in statistics}

Stepwise regression method

10 14.43 0.86 Variable X6 X7 X5 X8 X3 X4 X2 X9 X10 X11 P-value 0.11 0.13 0.15 0.23 0.33 0.48 0.50 0.51 0.63 0.88 9 14.50 0.86 Variable X7 X3 X8 X4 X6 X5 X10 X2 X9 P-value 0.00 0.00 0.00 0.00 0.04 0.14 0.37 0.51 0.52 8 14.04 0.86 Variable X4 X7 X8 X3 X6 X5 X10 X2 P-value 0.00 0.00 0.00 0.00 0.04 0.17 0.50 0.85 7 14.09 0.86 Variable X4 X7 X8 X3 X5 X6 X10 P-value 0.00 0.00 0.00 0.00 0.00 0.01 0.07 6 14.83 0.85 Variable X4 X7 X8 X3 X6 X5 P-value 0.00 0.00 0.00 0.00 0.00 0.00 5 15.66 0.82 Variable X4 X8 X7 X3 X6 P-value 0.00 0.00 0.00 0.00 0.03 4 17.86 0.80 Variable X3 X4 X8 X7 P-value 0.00 0.00 0.00 0.00 6 17.51 0.80 Variable X2 X3 X6 X7 X9 X10 P-value 0.00 0.00 0.00 0.00 0.00 0.08 GA 4 19.26 0.73 Variable X4 X7 X10 X6 P-value 0.00 0.01 0.01 0.98 Note: X2= SS, X3= HRT, X4= R%, X5= F/M, X6= SS× F/M, X7= SS× R%, X8= SS× HRT, X9= HRT× R%, X10= F/M× HRT, X11= F/M× R%.

(5)

Fig. 3. The results of evolutionary GM model.

Fig. 4. The results of GM model based on step-wise regression.

model can catch the whole tendency. In contrast, evolutionary GM model easily forecasts ﬂuctuations in the data series (Table 3). 3.2. Accuracy assessment

An artificial neural network (ANN) model is considered an effi-cient tool to predict the behaviors of a complex system with the characteristics of nonlinearity if sufficient observed data sets can be achieved. To prove the accuracy of the proposed model, a back-propagation neural (BPN) network was used to predict the COD concentration, then compared with the results of the proposed GM model.Fig. 5shows the architecture of the BPN model, which uses SS, HRT, R%, and F/M as input variables. The BPN network is trained with the generalized delta-learning rule, the gradient descent method to minimize error, sigmoid function as the acti-Table 3

Details of the BPN models used in this study.

Network architecture Optimization algorithm Optimum structure (input-hidden-output) Train cycles Random seed Learn rate Test period RMSE 4-3-1 3000 0.43 7 10 0.06114

RMSE means the root mean square error.

Fig. 5. The structure diagram of ANN used in this study.

Fig. 6. The results of artiﬁcial neural network based prediction model.

vation function, and the root mean square error (RMSE) to evaluate the performance of training and testing procedures.Fig. 6shows the result of a simpliﬁed BPN model. A comparison ofFigs. 3 and 6 reveals that the evolutionary GM model is a better performer in dynamic prediction than the simpliﬁed ANN model, based on the results of RMSE. The results also indicate the superiority of the proposed GM model because of their low data requirements. 3.3. Performance improvement

Decreasing the errors derived from the growth rate of pre-dicted sequences, this study not only adjusted the weights, but also optimized them by solving the nonlinear programming described above, with the number of variables restricted to four. Table 4 shows that the adjusted-R2_{varies from 0.71 to 0.82 as weights in}

Eq.(3)change, and the total errors can be reduced as the weights are determined by an optimization approach. Moreover, it appears that large residuals are associated with large monitoring values, Table 4

The results of different weight in calculating background coefﬁcient of GM model.

Weight Mean absolute

percentage error (MAPE) Adjusted-R2 p1 p2 p3 – 0.9 0.1 17.29 0.77 – 0.7 0.3 18.24 0.75 – 0.5 0.5 19.26 0.73 – 0.3 0.7 20.36 0.71 1/3 1/3 1/3 16.59 0.78 0.8 0.1 0.1 14.16 0.82 Z(1)_(k_{+ 1) = p} 1x(1)1 (k− 1) + p2x1(1)(k)+ p3x(1)1 (k+ 1).

(6)

Fig. 7. The results of evolutionary GM model with data transformed.

this study applied square root transformation to the variables. The modiﬁed GM model shows a discoverable improvement in accu-racy; e.g., the adjusted-R2_{increases to 0.85.}_{Figs. 7 and 8}_{show the} results of the GM(1, 4) prediction model with data transformed.

3.4. Sensitivity analysis of variables

To evaluate sensitivity response to the uncertainty of input vari-ables by using Monte Carlo simulation, R%, SS× F/M, SS × R%, and F/M× HRT were considered as the input variables. To quantify the confidence interval of the predicted values, 500 sample sizes were generated.Table 5presents the results of sensitivity analysis under different variable compositions. It shows the COD prediction is rel-atively sensitive to R% variation inferred by the narrowness of the confidence interval, but most sensitive to SS× F/M.Fig. 9portrays the confidence intervals for input items R% and SS× F/M.

Understanding the inﬂuences of observation and operational variables in prediction is useful for inexperienced operators while making a monitoring plan and deciding the quality requirements of monitoring variables. The conﬁdence intervals shown inFig. 6 provide an additional message for operators in applying the pro-posed model, together with error control. For example, operators are encouraged either to adopt more cautious actions if enormous variance results from these sensitive variables, or to increase reli-ability by obtaining monitoring data for special variables.

Fig. 8. The results of GM(1, 4) prediction model with data transformation.

Table 5

Sensitivity analysis under different variable compositions for GM(1, 4).

Variable compositions Conﬁdence interval of predicted COD concentration in efﬂuent Average upper boundary Average lower boundary R% 384.71 259.30 SS× R% 399.74 244.82 HRT, SS× F/M 483.58 164.05 F/M× HRT 516.85 131.99 SS× F/M 579.28 71.84 HRT, F/M× HRT 600.69 51.22 SS× R%, F/M × HRT 615.72 36.74 HRT, SS× F/M 663.12 −8.93 SS× F/M, SS × R% 678.15 −23.41 SS× F/M, F/M × HRT 795.26 −136.24 HRT, SS× F/M, SS × R%, F/M × HRT 977.96 −312.27

Fig. 9. Results of sensitivity analysis under different sampling compositions. (a) The influence of R% on confidence interval of predicted COD. (b) The influence of SS× F/M on confidence interval of predicted COD.

4. Conclusion

When a definite model is difficult to settle for a complex problem, a data-based prediction model becomes a considerable substitute. Genetic algorithm, capable of dealing with nonlinear optimization, is versatile at finding a behavior model based on a given database. Besides achieving accuracy predictions, this pro-posed approach in combining nonlinear programming with a GM model can simplify variable and parameter determination with lit-tle demand for operator experience in model building. Sensitivity

(7)

analysis results provide useful guides for operators to decide what future data should be collected for model reliability.

Acknowledgement

The authors express deep gratitude to the National Science Council in Taiwan for the ﬁnancial support (NSC 94-2211-E-242-004).

References

Akay D, Atak M. Grey prediction with rolling mechanism for electricity demand forecasting of Turkey. Energy 2007;32:1670–5.

Belanchea L, Valdes J, Comas JJ, Roda IR, Poch M. Prediction of the bulking phe-nomenon in wastewater treatment plants. Artiﬁcial Intelligence in Engineering 2000;1:307–17.

Belanche LA, Valdes JJ, Comas J, Roda IR, Poch M. Towards a model of input–output behaviour of wastewater treatment plants using soft computing techniques. Environmental Modelling & Software 1999;14:409–19.

Berthouex PM, Box GE. Time series models for forecasting wastewater treatment plant performance. Water Research 1996;30(8):1865–75.

Chang BR, Tsai HF. Forecast approach using neural network adaptation to support vector regression grey model and generalized auto-regressive conditional het-eroscedasticity. Expert Systems with Applications 2008;34:925–34.

Chen HW, Chang NB. Prediction of solid waste generation via grey fuzzy dynamic modeling. Resources Conservation and Recycling 2000;29:1–18.

Deng JL. Control problems of grey system. System Control Letters 1982;5:285–94. Gernaey KV, Van Loosdrecht MCM, Henze M, Lind M, Jørgensen SB. Activated sludge

wastewater treatment plant modelling and simulation: state of the art. Envi-ronmental Modelling & Software 2004;19:763–83.

Hanbay D, Turkoglu I, Demir Y. Prediction of wastewater treatment plant perfor-mance based on wavelet packet decomposition and neural networks. Expert Systems with Applications 2008;34:1038–43.

He RS, Hwang SF. Damage detection by a hybrid real-parameter genetic algorithm under the assistance of grey relation analysis. Engineering Applications of Arti-ﬁcial Intelligence 2007;20(7):980–92.

Hsu CC, Chen CY. Applications of improved grey prediction model for power demand forecasting. Energy Conversion and Management 2003;44:2241–9.

Lee DS, Vanrolleghem PA, Park JM. Parallel hybrid modeling methods for a full-scale coke wastewater treatment plant. Journal of Biotechnology 2005;115:317– 28.

Leephakpreeda T. Adaptive occupancy-based lighting control via grey prediction. Building and Environment 2005;40:881–6.

Lin ZC, Lin WS. Measurement point prediction of ﬂatness geometric tolerance by using grey theory. Precision Engineering 2001;25:171–84.

Marsili-Libelli S. Fuzzy prediction of the algal blooms in the Orbetello lagoon. Envi-ronmental Modelling and Software 2004;19(9):799–808.

Mao M, Chirwa EC. Application of grey model GM(1, 1) to vehicle fatality risk esti-mation. Technological Forecasting & Social Change 2006;73:588–605. Mjalli FS, Al-Asheh S, Alfadala HE. Use of artiﬁcial neural network black-box

mod-eling for the prediction of wastewater treatment plant performance. Journal of Environmental Management 2007;83:329–38.

Olson DL, Wu D. Simulation of fuzzy multiattribute models for grey relationships. European Journal of Operational Research 2006;175(1):111–20.

Sophie D, Adeline S. Estimation of parameters in incomplete data models deﬁned by dynamical systems. Journal of Statistical Planning and Inference 2007;137(9):2815–31.

Tarng YS, Juang SC, Chang CH. The use of grey-based Taguchi methods to determine submerged arc welding process parameters in hardfacing. Journal of Materials Processing Technology 2002;128:1–6.

Tenno R, Uronen P. State and parameter estimation for wastewater treat-ment processes using a stochastic model. Control Engineering Practice 1995;3(6):793–804.

Van Dongen G, Geuens L. Multivariate time series analysis for design and operation of a biological wastewater treatment plant. Water Research 1998;32(3):691–700. Wang CH. Predicting tourism demand using fuzzy time series and hybrid grey theory.

Tourism Management 2004;25(3):367–74.

Wang C, Hsu LC. Using genetic algorithms grey theory to forecast high technology industrial output. Applied Mathematics and Computation 2008;195(1):256– 63.

Wei YS. Efﬁcient and economical composting of sewage sludge for small and mid-scale municipal wastewater treatment plants. Ph.D. thesis. Research Center for Eco-Environmental Sciences, Chinese Academy of Sciences, Beijing, 2000. Wordena K, Wong CX, Parlitz U, Hornstein A, Engster D, Tjahjowidodo T, Al-Bender

F, Rizos DD, Fassois SD. Identiﬁcation of pre-sliding and sliding friction dynam-ics: grey box and black-box models. Mechanical Systems and Signal Processing 2007;21:514–34.

Zhou P, Ang BW, Poh KL. A trigonometric grey prediction approach to forecasting electricity demand. Energy 2006;31:2503–11.