Umer Khayam
email: umerkhayam@nwfpuet.edu.pk

Durre Nayab

Gul Muhammad Khan

S Ali Mahmud

Features Extraction of Growth Trend in Social websites using Non-Linear Genetic Programming

Keywords: Neuro Evolution, Artificial Neural Network, Cartesian Genetic Programming, Web Traffic Prediction, Web Portals, Future Demand

Nonlinear Cartesian Genetic Programming is explored for extraction of features in the growth curve of social web portals and establishment of a prediction model. Daily hit rates of web portals provide the measure of the growth and social establishment behavior over time. Non-linear Cartesian Genetic Programming approach also termed as CGPANN has unique ability of dealing with the nonlinear data as it provides the flexibility in feature selection, network architecture, topology and other necessary parameters selection to establish the desired prediction model. A number of socially established web portals are used to evaluate the performance of the model over a span of two years. Efficient performance is shown by the system keeping the fact in consideration that only single independent web portal data is used for training the network and the same network was used for the other web portals for their performance evaluation. The system performance is significantly good as the system selects only the desired features from the features presented as input and achieves an optimal network and topology that produce the best possible results.

INTRODUCTION

For the last few years a large number of web portals have been developed which intend to provide various kinds of valuable services to the internet users. Web portals such as social networks, online forums, job portals, blogging platforms, and news websites aim to provide uninterrupted free as well as paid services to their users. These web portals can also be used for business purposes or an advertising platform when its number of users increases and hit rate becomes high. The forecast of internet traffic is an important issue [1] and accurate forecasting of web portals helps in determining the current growth and future prospect of the portals. It also provides the business department an idea to identify the features and patterns [START_REF] Fausett | Fundamentals of Neural Networks: Architectures, Algorithms, and Applications[END_REF] that can lead to a higher hit rate in future, also an insight about the future demand and anomalies [START_REF] Krishnamurthy | Sketch-based change detection: methods, evaluation, and applications[END_REF] [START_REF] Jiang | Detecting network attacks in the internet via statistical network traffic normality prediction[END_REF] in the network can be detected. Software developers can get an idea about the types of browsers that are used to get access to the company website with the aid of forecasting. The business departments, developers and the network administrators do the task of forecasting internet traffic intuitively with the help of market information about the usual behavior and future number of visitors [START_REF] Papagiannaki | Long-term forecasting of Internet backbone traffic[END_REF]. This produces only a rough idea of what the future traffic will look like with a little day to day network administration. On the other hand, contributions from the areas of computational intelligence and artificial neural networks have replaced the intuition based forecasting methods [START_REF] Makridakis | Forecasting: Methods and Applications[END_REF] [START_REF] Hanke | Business Forecasting[END_REF]. The forecasting results obtained by these neural networks are accurate, reliable and less time consuming compared to the former methods.

Different methods have been used to forecast the internet traffic and several experiments have been performed on real-world datasets using these methods in order to determine their accuracy [1] [2]. In addition, different time scales have been used for forecasting purposes [2] and some methods perform well in short term forecasting, while others perform better in large term forecasts. By accurate prediction and forecasting, the services of web portals can be extended, diminished or altered in order to maintain a high number of visitors.

The main intent of this work is to present the application of artificial neural networks to web traffic forecasting for web portals.

LITERATURE REVIEW

TRAFFIC FORECASTING OF WEB PORTALS

Forecasting of web traffic is one of a critical task [START_REF] Papagiannaki | Long-term forecasting of Internet backbone traffic[END_REF] for researchers these days.

There are several applications of internet traffic forecasting that aims for efficient traffic engineering, anomaly detections and business tools development [1].

Various techniques have been used to predict the internet traffic accurately. Paulo Cortez used a neural network ensemble approach and two important adapted time series methods such as ARIMA and Holt-Winters [2] in order to determine their accuracy for predicting traffic in TCP/IP based networks. G. Rutka inspected the performance of traffic models for forecasting the future traffic variations as precisely as possible using multilayer perceptron and radial basis function networks based on measured traffic history [START_REF] Rutka | Neural network models for Internet traffic prediction[END_REF]. Machine learning methods have been used to analyze web traffic using decision trees [START_REF] Piramuthu | On learning to predict web traffic[END_REF]. A multi scale decomposition approach was used for real time traffic prediction that outperforms traffic prediction using neural network approach [START_REF] Mao | Real-time network traffic prediction based on a multiscale decomposition[END_REF] and gives comparatively better results. Many internet service providers use multi-protocol label switching to establish fill mesh of MPLS between pairs of the routers in a network in order to optimize the bandwidth resources in the network [START_REF] Davie | MPLS: technology and applications[END_REF]. Even if MPLS is not used, with the knowledge of future demand of traffic matrix, the traditional allocation of the routing protocol weights can be done more efficiently [1]. The forecasted web traffic cab be used for anomaly detection [11] [12] in the network by comparing the actual traffic to the forecasted traffic in the network. The forecasted web traffic can be used by the business departments [1] of the company to get an idea about the current popularity and future demand of their web services and applications. The statistical data about the visitors on a particular web portal contain patterns that are viable in determining the future demand. By analyzing statistical data, the factors effecting hit rate of the web portals can be identified, an insight about the future hit rate, types of visitors, and survival of the web portals can be predicted. Several factors have been identified in past that effect the number of visitors on these portals [1]. These factors include accessibility, browsers compatibility, trust level, flexibility, productive and organized content, and the targeted age, group and gender etc [2]. The business department and the developers of the company can take decisions based on the forecasted results in order to get increased traffic in the future.

Cartesian genetic programming and neural network models are used for forecasting purposes. Today the numbers of structures of the artificial neural networks and training algorithms that are applied in their learning process have evolved and these models are also used to forecast future network traffic. Similarly these models can be used to forecast the traffic patterns of web portals.

CARTESIAN GENETIC PROGRAMMING (CGP)

The idea of Cartesian Genetic Programming (CGP) was proposed in 1999 by Julian F Miller [START_REF] Rothermich | Studying the emergence of multicellularity with cartesian genetic programming in artificial life[END_REF], that is a genetic programming approach with a two dimensional graphical arrangement of its functional units i.e. nodes. These programs are represented as directed acyclic structures operating in a feed forward fashion. CGP has two formats of network representation i.e. phenotype and genotype. The physical format of array of interconnected nodes represents the phenotypic structural space and its genotype can be represented as an array of integers i.e. genes. The nodes are the functional units that exhibit a node function, inputs, and outputs [START_REF] Miller | Cartesian Genetic Programming[END_REF]. CGP provide a general platform for evolving the hybrid structure of any number of networks in any order. It provides a complete Cartesian architecture for their interconnectivity patterns producing less hybrid structures from a pole of networks. CGP is explored in a range of applications producing interesting results [START_REF] Miller | Cartesian Genetic Programming[END_REF].

CARTESIAN GENETIC PROGRAMMING EVOLVED ARTIFICIAL NEURAL NETWORK

The strategy employed for the evolution of an ANN plays a vital role in ANNs and hence is a major concern for the researchers these days. The evolutionary strategy used for the evolution of ANN model proposed in this research is CGP. CGPANN is signal processing system that is based on some of the known organizing principles of the human brain [START_REF] Fausett | Fundamentals of neural networks[END_REF] [16] [START_REF] Zurada | Introduction to artificial neural systems[END_REF]. CGP evolve the ANN with its unique architecture that makes it computationally cost effective and efficient. They are computational models that are capable of machine learning and pattern recognition. These systems are represented by a number of independent, simple processors called neurons which are interconnected with each other by weighted connections [START_REF] Ding | Neural network based models for forecasting[END_REF]. Figure 1 shows a typical CGPANN genotype and phenotype with inputs (I0, I1, I2, ... I9), outputs (O0, O1, O2, ... O9) , active nodes (0, 1, 2, 4, 5, 7, 8 & 9), inactive nodes (3 & 6) and weights (w0, w1, w2……., w17). The arity (number of inputs per node) of the network is 2 and the number of inputs to the system is 10. The string of numbers underneath the figure is the genotype arranged in boxes represent individual node in terms of its function, inputs and the weights associated with the inputs respectively. The box having dotted lines represents the inactive nodes. The generic expressions for the CGPANN model can be given in the following equations. Equation (1) shows the system outputs as the summation of system inputs ( ) when the weights and node functions are not associated with them. Equation ( 2) and (3) shows the network outputs with the weights (wi) and functions ( ) are associated with the system inputs.

∑

(1)

∑ (2) (∑ ) (3)
Equation ( 4) and [START_REF] Nayab | Prediction of Foreign Currency Exchange Rates Using CGPANN[END_REF] shows the range (i) of the inputs, functions and weights respectively. As can be seen from Eq (4) the range (i) belongs to natural numbers (N) where N begins from 1 to NT. The weights belong to real numbers and are chosen in the range [-1, 1].

(5)

Equation ( 6) represent the single generation obtained from the network in terms of inputs (I), system outputs (yi's) and ultimate output (Op).

(6)

Equation ( 7) shows the ultimate system output in terms of individual system outputs. Equation (8) shows the individual system output in terms of nodal outputs (yi) and system inputs (I's) with weights (wi's) associated with them.

∑ (8) Oi= ∑ (9)

MUTATION IN CGPANN

The optimal network of the CGPANN is achieved during the process of evolution. During evolution mutation in the genotype takes place and optimal network is achieved by the selection and promotion of the fittest genotype. Figure 2 (a,b,c) show the mutation process of the CGPANN network during evolution process.

Figure 2a shows the original phenotype and genotype of the network and Figure 2b and 2c demonstrate the process of mutation in the function (f2) and input (I3) genes respectively.

THE FORECASTER MODEL

The learning process of our forecaster model consists of adaptive modifying of the connection weights to improve the overall functionality of the neural network as the parallel signal processing system. The system learns the patterns and selects the optimal features in the historical data to predict the future values.

The performance of the forecasting model is measured by Mean Absolute Percentage Error (MAPE) as it is a common metric in forecasting applications [START_REF] Papagiannaki | Long-term forecasting of Internet backbone traffic[END_REF].

The forecaster model consists of single row and multiple columns of neurons trained by altering the values of the connections between the neurons. A nonlinear CGP is used for translating the network as it handles the nonlinear data such as that of web portals efficiently. The statistical data of six web portals for a period of one year is fed to the forecasting model and model is trained. The network is capable of feature selection i-e selecting optimal (not all) number of inputs and nodes for evolving the final network. The neural network model is trained on a one year hourly spaced datasets of Linkedin.com. The trained model is then used to forecast the internet traffic of web portals other web portals. The mean absolute percentage error is calculated between the forecasted values and the actual values.

EXPERIMENTAL SETUP

A collection of time ordered dataset is taken for six web portals including Linkedin.com, Time.com, Tumblr.com, Answers.com, Hubpages.com and Collegehumor.com for a period of two years starting from October 2011 to October 2013. The data is obtained and verified by Quantcast.com that is a web portal used for monitoring the visitor hit rate on the portals available on the internet. The obtained data was the average daily hit rates of the six specified web portals for a period of two years.

During the training phase daily global hit rates of LinkedIn.com is used for training the system. Initially a random population of ten networks is produced for five independent seeds. The system parameters are defined and initialized. The number of offspring per generation is 9 under the 1+λ evolutionary strategy where λ is set to be 9. The initially generated genotype is mutated to produce nine offspring. The best genotype among these genotypes is selected and mutated again. The process is repeated unless the best network is achieved. The mutation rate (μr) is set to be 10%. The number of system inputs is set to be ten and that is the daily hit rates data and the arity of the system i.e. number of nodal inputs is set to be 5. Log sigmoid is taken as the activation function. These system parameters are chosen based on previous performance of CGPANN [START_REF] Nayab | Prediction of Foreign Currency Exchange Rates Using CGPANN[END_REF][START_REF] Pacelli | An Artificial Neural Network Model to Forecast Exchange Rates[END_REF] and evolutionary performance. The network is trained for variable number of inputs and outputs. The network takes seven days and fourteen days instances of the daily dataset as inputs to the network and forecasts the hit rate of the web portal for the next single day, seven days and fourteen days as outputs of the network.

Performance of the system is evaluated with MAPE (Mean Absolute Percentage Error) that is calculated by comparing the actual value of the hit rate with that estimated by the model. The mathematical expression for MAPE is given below:

∑ ( )

Where LAi is the actual value, LFi is the estimated value and N is the number of days.

RESULTS AND ANALYSIS

During the training phase the system is been trained for one million generations and the best trained network is then tested for its performance. The system is trained on one web portal data i.e. LinkedIn.com and is tested on a different dataset of five web portals. The data with known hit rate values are estimated during the testing phase and compared with the actual values to evaluate the performance of the model. The CGPANN forecaster model proposed in this work has been trained using daily averaged web traffic data of LinkedIn.com for one year, starting from October 2011 to October 2012. The results for the training session of each network are given in the following Table I. The output of the system is evaluated from the average of five independent evolutionary runs for each network with mentioned combination of inputs and outputs. The fittest network is achieved for the network with 7 inputs and 1 output during the training phase with the MAPE value of 1.1149. the mentioned combination of inputs and outputs. The proposed forecaster model takes seven days hit rates as input to the system and predicts the eighth day hit rate in the first case. In the second case it takes seven daily hit rates as input and predicts seven future daily hit rates as output of the system. In all the combinations of the inputs and outputs the forecaster model showed higher accuracy hence the model is proficient both in terms of single instant prediction and multiple instants predictions.

CONCLUSION AND FUTURE WORK

We have explored nonlinear CGP for the implementation of forecasting model for global web portals growth analysis and extracting the prominent features. A number of socially established web portals are analyzed for their growth process. The performance of the system revealed that the system is robust that learns the trends and extracts the optimal features responsible for the growth rate of these portals. The network is capable of feature selection i-e selecting optimal (not all) number of inputs and nodes for the evolution of the final network. The proposed system has the ability to obtain an optimal set of features, number of nodes and connections paradigm and morphology for the best possible prediction model for the task at hand. Further work can explore the social behavior on individual portals including: probing the posts and updates, advertisement on web portals for improving the news feeds, analyzing the intent and emotions in user updates, company survival capabilities analysis and business success or failure of start-up firm analysis.

9 FIGURE 1 :

 91 FIGURE 1: A TYPICAL CGPANN PHENOTYPE AND GENOTYPE

FIGURE 2A :FIGURE 2B :FIGURE 2C :

 2A2B2C FIGURE 2A: ORIGINAL CGPANN PHENOTYPE Genotype1= f0 I 0 w 0 I 1 w 1 I 2 w 2 , f1 f 0 w 7 I 4 w 4 I 2 w 6 , f2 f 0 w 7 I 4 w 8 I 3 w 9 , f3 f 0 w 10 f 1 w 11 I 0 w 12 I0

Table 1 :

 1 Training results of CGPANN for web portal forecasting model

	Nodes	7 in 1 out	7 In 7 out	14 In 7 out
	50	0.11675	0.13302	0.132575616
	100	0.11686	0.13457	0.132009816
	150	0.11493	0.13291	0.13219517
	200	0.11581	0.13402	0.132067988
	250	0.11624	0.13377	0.133859813
	300	0.11557	0.13290	0.132211834
	350	0.11666	0.13318	0.132714417
	400	0.11619	0.13349	0.133092481
	450	0.11925	0.13347	0.13273484
	500	0.11659	0.13286	0.13262616
	For evaluation and testing, a new dataset is fed to the network for its validation.
	The new dataset is the statistical data of the visitor hit rate of Time.com, Tum-
	blr.com, Answers.com, Hubpages.com and Collegehumor.com for a period of two
	years from ranging from October 2011 to October 2013. The testing phase re-
	sults for the hit rates forecasting model are tabulated in Table II, III and IV for

  The best MAPE values achieved for each web portal are highlighted in each table respectively. The best MAPE value achieved by the model is 9.053% given in Table2where as other works done the web portal forecasting has achieved accuracies in the ranges such as 12-23% [1], 13-22%[2], 72.74 and 12.04 %[START_REF] Gluszek | Web traffic prediction with artificial neural networks[END_REF].

Table 2 :

 2 Testing results for one day hit rate as input with seven day's data as output

	No. Of Nodes Times	Tumblr	Answers	Hubpages	Collegehumor
	50	0.1463	0.09053	0.14202	0.09647	0.1553
	100	0.1433	0.09261	0.14445	0.09671	0.1541
	150	0.1549	0.09426	0.14266	0.10058	0.1692
	200	0.1537	0.09373	0.14592	0.10016	0.1652
	250	0.1603	0.09414	0.14548	0.10192	0.1730
	300	0.1430	0.09072	0.14254	0.09622	0.1507
	350	0.1469	0.09361	0.14625	0.09828	0.1565
	400	0.1545	0.09301	0.14100	0.09983	0.1677
	450	0.1686	0.10202	0.15097	0.10875	0.1873
	500	0.1424	0.09478	0.14661	0.09819	0.1542

Table 3 :

 3 Testing results for seven days hit rates as input with seven days data as output

	No. Of Nodes Time	Tumblr Answers	Hub pages	College humor
	50	0.1442	0.1159	0.19272	0.10967	0.15613
	100	0.1426	0.1152	0.19391	0.10884	0.15269
	150	0.1270	0.1141	0.18955	0.10495	0.13731
	200	0.1496	0.1191	0.19283	0.11337	0.16625
	250	0.1571	0.1201	0.19474	0.11554	0.17485
	300	0.1447	0.1165	0.19327	0.11035	0.15741
	350	0.1530	0.1186	0.19458	0.11311	0.16570
	400	0.1351	0.1151	0.19247	0.10721	0.14410
	450	0.1112	0.1141	0.19209	0.10121	0.11393
	500	0.1121	0.1134	0.19011	0.10173	0.11769

Table 4 :

 4 Testing results for 14 days hit rates as input with 7 days data as output

	No. Of Nodes Time	Tumblr	Answers	Hub pages	College humor
	50	0.12795	0.11923	0.19835	0.10806	0.14111
	100	0.13847	0.11873	0.19897	0.10981	0.15333
	150	0.15067	0.12162	0.20324	0.11404	0.16946
	200	0.13021	0.11689	0.19627	0.10698	0.14178
	250	0.09451	0.11605	0.19751	0.09905	0.09592
	300	0.15046	0.12001	0.20089	0.11343	0.16898
	350	0.14978	0.12036	0.20373	0.11343	0.16571
	400	0.14785	0.12161	0.20692	0.11215	0.16292
	450	0.15392	0.12095	0.20030	0.11530	0.17346
	500	0.14222	0.11855	0.19732	0.11093	0.15814


Petros-Fotios Alvanitopoulos

Ioannis Andreadis

Nikolaos Georgoulas

Michalis Zervakis
email: michalis@display.tuc.gr

Nikolaos Nikolaidis
email: niko-laos.nikolaidis@enveng.tuc.gr

Solar radiation time-series prediction based on Empirical Mode Decomposition and Artificial Neural Networks

This paper presents a new model for daily solar radiation prediction. In order to capture the hidden knowledge of existing data, a time-frequency analysis on past measurements of the solar energy density is carried out. The Hilbert-Huang transform (HHT) is employed for the representation of the daily solar irradiance time series. A set of physical measurements and simulated signals are selected for the time series analysis. The empirical mode decomposition is applied and the adaptive basis of each raw signal is extracted. The decomposed narrowband amplitude and frequency modulated signals are modelled by using dynamic artificial neural networks (ANNs). Nonlinear autoregressive networks are trained with the average daily solar irradiance as exogenous (independent) input. The instantaneous value of solar radiation density is estimated based on previous values of the time series and previous values of the independent input. The results are promising and they reveal that the proposed system can be incorporated in intelligent systems for better load management in photovoltaic systems.

Introduction

Today there is an obvious need for eco-friendly and sustainable energy sources. It is well known that in developed world the daily life depends on the continuous supply of energy. Unfortunately, the exponential increase in energy demand compromises the ability of future generations to meet their needs. The provision of fossil fuel cannot be taken for granted. In recent years many European govern-ments introduced subsidy programs in order to encourage the utilization of alternative energy resources with high efficiency and environmentally sound.

Nowadays, there are many options for renewable energy sources. One of the most attractive and promising resource is the solar irradiance. Solar energy usage is small globally, but it is commonly accepted that it will increase substantially in the near future.

Solar energy is one of the most considerable design parameters in many research fields such as bioclimatic design of buildings, analysis of the performance of photovoltaic systems, agriculture science and meteorology. Moreover, in the case of autonomous energy systems with solar power source an automatic and efficient management of solar energy is required in order to minimize the strain of storage components and increase their lifetime.

Due to the widely use of solar energy there is an obvious need for a more precise modeling and close prediction of the solar irradiance future trend. Several studies for solar prediction exist in the literature. These studies involve empirical models [START_REF] Batlles | Empirical modelling of hourly direct irradiance by means of hourly global irradiance[END_REF][START_REF] Loutzenhier | Empirical validation of models to compute solar irradiance on inclined surfaces for building energy simulation[END_REF], statistical approaches [START_REF] Safi | Prediction of global daily solar radiation using higher order statistics[END_REF][START_REF] Kaplanis | A model to predict expected mean and stochastic hourly global solar radiation I(h;nj) values[END_REF] and intelligent system approaches [START_REF] Moustrisa | Use of neural networks for the creation of hourly global and diffuse solar irradiance data at representative locations in Greece[END_REF][START_REF] Rehman | Artificial neural network estimation of global solar radiation using air temperature and relative humidity[END_REF].

This paper presents a new approach for solar radiation prediction by taking into account the past measurements of the series. The HHT operates on the original signal in order to extract the different spectral components. Afterwards, nonlinear autoregressive neural networks are trained with the average daily solar irradiance as exogenous input. The instantaneous values of solar radiation density are predicted for each of the partial component and are summarized in order to estimate the original signal's prediction.

The paper is organized as follows: Section 2 introduces the HHT. Section 3 refers to the nonlinear autoregressive networks, while section 4 deals with the proposed model implementation. Section 5 presents the experimental results.

The Hilbert Huang Transform

The HHT is an appropriate tool for analyzing non-stationary and nonlinear data. The HHT, a NASA's designated name, is proposed by Huang et al [START_REF] Huang | The empirical mode decomposition and the Hilbert spectrum for nonlinear and non-stationary time series analysis[END_REF]. The essential difference between HHT and most currently data decomposition methods is that HHT provides an adaptive analysis of data. The term adaptivity refers to the requirement for an adaptive basis. There is no a priori reason to believe that an arbitrarily selected basis is able to describe the underlying processes of any given phenomenon. This algorithm comprises Empirical Mode Decomposition (EMD) and Hilbert Spectral Analysis (HSA). The combination of EMD and HSA provides a consistent and physically meaningful definition for instantaneous frequency and amplitude. First, the EMD process decomposes any complicated set of data into a finite and often small number of intrinsic mode function components. This decomposition is adaptive since it is based on local characteristics of the original signal and it can be applied to nonlinear and nonstationary signals. Second, the HSA is applied to the generated components in order to provide a meaningful time-frequency-energy distribution of time series.

Empirical Mode Decomposition (EMD)

The EMD step is based on the assumption that any time series consists of different, intrinsic modes of oscillation. Each mode of oscillation is derived from the data and fulfills the following conditions:

1 In the whole data set, the number of extrema and the number of zerocrossings must be equal or differ at most by one. 2 At any point, the mean value of the envelope defined by the local maxima and the envelope defined by the local minima is zero. An intrinsic mode function (IMF) represents a simple oscillatory mode similar to a component in the Fourier-based simple harmonic function. IMFs constitute an a posteriori adaptive basis. The decomposition of any time series is described below.

First, all the local extrema of the series x(t) are extracted. The local maxima are connected by a cubic spline to produce the upper envelope u(t). Similarly the lower envelope v(t) is created with the connection of the local minima. The mean of these two envelopes is designated as m1, and the difference between x(t) and m1 is the first component h1.

(

(2) Usually, h1 does not satisfy all the requirements of an IMF. Therefore, the previous process, entitled "sifting process", is repeated more times. In the subsequent sifting process, h1 is treated as the data, then

(3) where m11 is the mean of the upper and lower envelopes of h1. After repeated sifting, up to k times which is usually less than 10, h1k, given by ( 4) is designated as the first IMF component c1 from the data, or (5) Typically, c1 will contain the shortest-period component of the signal. The c1 component is removed from the original signal in order to obtain the residue (6) The residue r1, which contains longer-period components, is treated as the new data and subjected to the same sifting process as it has been described above. This procedure can be repeated to obtain all the subsequent rj functions as follows:

(7) The sifting processs is terminated by either of the following predetermined criteria: either the component cn or the residue rn becomes so small that it is less than a predetermined value of consequence, or the residue, rn, becomes a monotonic function, from which no more IMF can be extracted. Using the above equations the original data is the sum of the IMF components plus the final residue (8)

Hilbert Spectrum Analysis (HSA)

In HSA the Hilbert transform (HT) is applied in each extracted IMF to obtain the representation of the original signal in time-frequency-energy domain. The application of HT in each IMF yields sharp frequency and time localization. For an arbitrary function, x(t) its Hilbert transform, y(t), is defined as [START_REF] Gao | NARMAX time series model prediction: feedforward and recurrent fuzzy neural network approaches[END_REF] and the analytic signal z(t) is described by (10) where [START_REF] Moirogiorgou | High Frequency Monitoring System for Integrated Water Resources Management of Rivers[END_REF] the instantaneous frequency is given by

f ( 12
)
and the original signal can be expressed by ( 13) where (13) represents the general form of Fourier equation ( 14)

3 Nonlinear Autoregressive Network with Exogenous Inputs

NARX Networks

Several processes in real world are described by time series. The ANNs have been widely used in the scientific field of time series prediction due to their inherent nonlinearity and high robustness in noise. It is well-known that ANN represents a nonparametric approach since there is no need to have any knowledge of the underlying processes that generate the original signal.

Typically, the challenge task of time series prediction can be expressed as finding the appropriate function I so as to acquire an estimate of the time series y at time t+D (D=1,2 …) given the past values of y up to time t, plus the values of exogenous input x:

(15) where y(t) and x(t) represent the values of y and x in time t respectively. The variables and are the lags parameters of model and in case of D>1 we have the multi step ahead prediction of time series y.

NARX networks are recurrent neural architectures and it has been demonstrated that they are powerful in theory and at least equivalent to Turing machines. One can utilize those models rather than conventional recurrent networks without any computational loss [START_REF] Siegelmann | Computational capabilities of recurrent narx neural networks[END_REF]. The main advantages of NARX networks are: more effective learning of models and faster convergence than in other ANNs [START_REF] Gao | NARMAX time series model prediction: feedforward and recurrent fuzzy neural network approaches[END_REF].

Proposed Method

Solar radiation dataset

In this paper, a study of forecasting half daily solar radiation is presented by taking into consideration past values of the solar radiation time series. The dataset contains solar radiation measurements provided by the Joint Research Centre (JRC) Renewable Energies unit of the European Commission [10]. Cubic spline interpolation is utilized in order to increase the size of data. Thus, from a 15minute sample frequency is created a 1-minute sample frequency time series. Each daily artificial solar irradiance signal is a vector of length 840 which represents the daily 1-minute time interval solar irradiance samples. The geographic parameters latitude and longitude were set to 35°27'29" and 24°9'5", respectively. These parameters refer to Koiliaris river region in Crete, where the solar radiation prediction is studied to manage the energy demands of the CyberSensors monitoring system [START_REF] Moirogiorgou | High Frequency Monitoring System for Integrated Water Resources Management of Rivers[END_REF]. The daily solar power density samples, which are utilized in our analysis, are records that refer to the most adverse solar radiation conditions. This data is taken in December, which is the month that provides the lowest amount of solar energy.

Solar radiation analysis

In the present study the "clear-sky" values and the data for "real" conditions were taken into account. The former refers to values of irradiance where the sky is completely free of clouds. This signal is valuable since it provides the maximum values of solar irradiance that may appear at a given time in the examined month. The latter refers to values of solar irradiance with average cloud cover for the given month. In Fig. 1 the "clear-sky" and the time series with average values of solar irradiance in December are presented for the examined region. A typical solar radiation power density for a cloudy day is shown in Fig. 2.

Fig. 2. Solar radiation power density for a cloudy day in December.

The decomposition of each examined signal is implemented by means of EMD. An example of the EMD is presented in the following figure, where the solar power density signal from the previous figure is analyzed into a set of IMFs. It should be mentioned that the generated basis is an a posteriori basis created based on the local characteristics (local maxima and minima) of the time series. From Fig. 3 it can be noticed that the initial signal consist of five IMFs and a residue.

The NARX models

Nonlinear autoregressive with exogenous inputs networks were used for the solar radiation prediction. For each extracted intrinsic mode function a NARX network is trained and used for the prediction of the specific mode of oscillation. The NARX network is a two-layer feed forward artificial neural network. In the hidden layer a sigmoid transfer function is utilized and a linear transfer function is set to the output layer. In order to predict the future values of the time series y(t), which in our case represents the solar radiation intrinsic mode function, the past values of time series x(t) are also required. In the present study the clear-sky solar irradiance power density is treated as the exogenous input x(t). Fig. 4 shows the topology of the NARX network, where the values dy and dx are the delays of y(t) and x(t) respectively.

Results

NARXs models are trained taking into account the solar radiation time series, which is derived from measurements of twenty five days. Five days records are utilized to assess the performance of the algorithm. The training function is set to Levenberg-Marquardt algorithm and the performance metric is the Mean Squared Error (MSE). The maximum horizon, which is taken, is 420 time steps ahead. The performance is estimated by taking into consideration the root-mean squared error (RMSE) and the maximum error. Table 1 shows the simulation results. The number of hidden neurons is set to 10 and the number of delays is 30. The validation and test data points are set to 15% of the original data, while the rest 70% is used for training.

Fig. 5 shows the simulation results of solar radiation for a half day horizon. As it can be seen from the below figures, solar radiation prediction is a complex task due to the extremely uncorrelated high frequency components. These high frequency components represent clouds and atmosphere attenuation, which are extremely difficult to predict using only past information. Thus, additional inputs are required in order to optimize the simulation results.

Conclusions

In this paper the HHT is used for the decomposition of solar irradiance time series in order to predict solar radiation within a time horizon of half day ahead. From the results it can be seen that a more accurate prediction is required for the higher frequency components of the solar radiation time series. Taking into account extra data such as images from satellites, additional information regarding the impact of cloudiness on ground irradiance will be obtained. Moreover, deferencing of solar irradiance time series may generate stationary signals which are relatively easy to predict. In the nearby future, research will continue in this line by following the aforementioned approaches with ultimate goal to integrate this technique in an intelligent Energy Consumption Controller.

Fig. 1 .

 1 Fig. 1. Solar irradiance plot in December for the examined region with 25 o angle.

Fig. 3 .

 3 Fig. 3. Decomposition of a solar radiation signal into IMF components.

Fig. 4 .

 4 Fig. 4. Topology of the NARX network.

  Fig 6 shows an alternative view of the difference between the desired and predicted solar radiation values for the first day.

Fig. 5 .

 5 Fig. 5. Simulation of solar radiation for a half day maximum horizon.

Fig. 6 .

 6 Fig. 6. Desired and predicted solar radiation values for the first day.

Table 1 . RMSE from simulation results.

 1

		1st day	2nd day	3rd day	4th day	5th day
	RMSE (W/m 2 )	26.891	21.732	31,009	39,658	35,729
	MAX Error	97.847	60.130	120.510	129.815	128.858

Acknowledgements

This work is elaborated through an on-going THALES project (CYBERSENSORS -High Frequency Monitoring System for Integrated Water Resources Management of Rivers). The Project has been co-financed by the European Union (European Social Fund -ESF) and Greek national funds through the Operational Program "Education and Lifelong Learning" of the National Strategic Reference Framework (NSRF) -Research Funding Program: Thales. Investing in knowledge society through the European Social Fund.


Mehreen Rehman
email: mehreen@nwfpuet.edu.pk

Jawad Ali
email: jawad.ali@nwfpuet.edu.pk

Gul Muhammad Khan

Sahibzada Ali Mahmud
email: sahibzada.mahmud@nwfpuet.edu.pk

Extracting Trends Ensembles in Solar Irradiance for Green Energy Generation using Neuro-Evolution

Keywords: Solar irradiance Forecasting, Neural Networks, Cartesian Genetic Programming (CGP), Neuro-evolution, Time Series Prediction

Globally, there are variations in climate, there is fossil fuel depletion, rising fossil fuel prices, increasing concern regarding energy security, and awareness about the environmental impacts of burning fossil fuels. These factors lead to a growing interest around the world in green and renewable energy resources, solar energy being a common one. A Neuro-evolutionary approach is explored to extract the trend ensembles in the solar irradiance patterns for renewable electric power generation, using the data taken from stations in Al-Ahsa, Kingdom of Saudi Arabia. The algorithm, based on Cartesian Genetic Programming Evolved Artificial Neural Network (CGPANN) was developed and trained for hourly and 24-hourly prediction, using the solar irradiance value as the input parameter. It was tested to predict solar irradiance on hourly, daily, and weekly basis. The proposed technique is 95.48% accurate in solar irradiance prediction.

Introduction

Most places on the earths surface receive clean and abundant solar energy, free of cost, throughout the year. Owing to the rising fossil fuel costs and the degradation of atmosphere by these fossil fuels, there is a dire need for economical and effective harnessing and utilization of solar energy. Because the world has to pay the cost of losing its energy resources at an alarmingly fast pace, there is an increasing reliance of power generation on renewables, such as solar energy. Worldwide efforts are being taken to increase the capacity of renewable power production. However, the intermittent nature of renewables makes it a challenging task to integrate them into any system. Accurate solar irradiance is therefore essential in improving the efficiencies of solar and wind energy based applications, obtaining optimized power production [START_REF] Koeppel | Using storage devices for compensating uncertainties caused by non-dispatchable generators[END_REF], overcoming the challenges of integrating them in the power grid and for the efficient management and operation of solar thermal energy plants. Solar radiation data gives us information about the amount of energy that hits the earths surface during a particular time frame and is required for research on effective solar energy utilization. Because this data is not readily available, and because the intensity and availability of solar radiation is dependent upon several other environmental factors, alternate ways of generating the data must be developed [START_REF] Dorvlo | Solar radiation estimation using artificial neural networks[END_REF].

Artificial Neural Networks is an emerging technology for handling complex, practical problems, such as forecasting. Neural networks are capable of making computationally efficient, faster and more practical predictions than any of the conventional methods.

Our work mainly focuses on introducing a new and efficient solution based on the Neuro-evolutionary technique termed as the Cartesian Genetic Programming Evolved Artificial Neural Network (CGPANN) [START_REF] Khan | Efficient representation of recurrent neural networks for markovian/non-markovian non-linear control problems[END_REF] to extract the trend ensembles in the solar irradiance patterns for renewable power generation. Neuro-Evolution is the phenomenon that involves the artificial evolution of the entities of an Artificial Neural Network (ANN). In CGPANN, Neuro-Evolution exploits the powerful structural properties of Cartesian Genetic Programming [START_REF] Miller | Cartesian genetic programming[END_REF] and the functional properties of an ANN. This leads to the evolution of all the network parameters; node weights, system inputs, node inputs, node functions, network topology as well as the system outputs [START_REF] Khan | Efficient representation of recurrent neural networks for markovian/non-markovian non-linear control problems[END_REF]. The research solution discussed here is specific only to CGPANN. The proposed technique uses past solar irradiance data to estimate the hourly, 12-, 24-, 48-and 168-hourly solar radiation patterns.

Literature Review

Computational Intelligence Techniques for Solar Irradiance Forecasting

Estimation methods used in the literature for forecasting solar radiation include stochastic, analytic [START_REF] Dagestad | Estimating global radiation at ground level from satellite images[END_REF], and artificial neural networks [START_REF] Mohandes | Use of radial basis functions for estimation monthly mean daily solar radiation[END_REF].The analysis by Mubiru and Banda [START_REF] Mubiru | Estimation of monthly average daily global solar irradiation using artificial neural networks[END_REF] shows that owing to the nonlinear, non-stationary nature of solar radiation; ANNs are thought to be superior to empirical methods in estimating the global solar radiation. Autoregressive (AR) and autoregressive moving average (ARMA) are models frequently used for the estimation of solar irradiance [START_REF] Reikard | Predicting solar radiation at high resolutions: a comparison of time series forecasts[END_REF][START_REF] Torres | Forecast of hourly average wind speed with ARMA models in Navarre (Spain)[END_REF]. Angstrom [START_REF] Angstrom | On the computation of global solar radiation from records of sunshine[END_REF] was the pioneer to apply Regression models based on sunshine duration. Ahhi et al. [START_REF] Ahhi | Prediction of monthly average daily global radiation in Al Ain city-UAE using artificial neural networks[END_REF] have used various input parameters of an ANN to predict solar radiation.

Kalogirou Model applies the standard back propagation learning algorithm. The model uses the hourly records for one year to predict hourly solar radiation. Other ANN models are; The Mohandes Model [START_REF] Mohandes | Use of radial basis functions for estimation monthly mean daily solar radiation[END_REF]; The Kemmoku Model [START_REF] Koeppel | Using storage devices for compensating uncertainties caused by non-dispatchable generators[END_REF]; The Reddy Model [START_REF] Sozen | Solar-energy potential in Turkey[END_REF]; The Sozen Models [START_REF] Cao | Forecast of solar irradiance using recurrent neural networks combined with wave analysis[END_REF] that has explored Scaled conjugate gradient (SCG), Pola-Ribiere conjugate gradient (CGP), and Levenberg-Marquardt (LM) learning algorithms and logistic sigmoid transfer function; The Cao Model [START_REF] Soares | Model-ling hourly diffuse solar radiation in the city of Sao Paulo using artificial neural network technique[END_REF] that establishes a recurrent BP network for forecasting of solar irradiance; The Soares Model [START_REF] Linares-Rodriguez | An artificial neural network ensemble model for estimating global solar radiation from Meteosat satellite imag-es[END_REF] that develops a neural network exploring feature determination and pattern selection techniques. Linares-Rodriguez et. al. [START_REF] Bacher | Online short-term solar power forecasting[END_REF] explores an articial neural network ensemble model for estimating global solar radiation from Meteosat satellite images. This model uses clear-sky estimates and satellite images as input variables. Tymvios et al., Angela et. Al [START_REF] Hocaoglu | Hourly solar radiation forecasting using optimal coefficient 2-D linear filters and feed-forward neural networks[END_REF] have introduced a single-parameter model for solar irradiance prediction. Lorenz et al., [START_REF] Cao | Study of hourly and daily solar irradiation forecast using diagonal recurrent wavelet neural networks[END_REF] use numerical weather predictions (NWPs) as input parameters to NNs to estimate global irradiance. Hocaoglu et al., [START_REF] Jain | Prediction of Global Solar Radiation Using Artificial Neural Network[END_REF] study and compare feed-forward NNs with seasonal AR models. Cao and Lin, [START_REF] Azadeh | An integrated artificial neural networks approach for predicting global radiation[END_REF] combined NNs with wavelets to predict next day hourly values of global irradiance, incorporating various meteorological parameters as inputs to the models. Jain and Goel, [START_REF] Yang | Improving estimation of hourly, daily, and monthly solar radia-tion by importing global data sets[END_REF] developed a multilayer feed forward neural network model for solar radiation forecasting in India. For the input to the network, the geographical, solar and meteorological parameters were taken from different areas of India. Azadeh et. al. [25] proposed an ANN model based on multilayer perceptron approach, incorporating all climatological and meteorological factors as input variables for prediction of solar global radiation. The results of the integrated ANN-MLP model have been compared to those of angstroms model [START_REF] Angstrom | On the computation of global solar radiation from records of sunshine[END_REF]. There is a notable decrease in the MAPE value from 14.78% by the Angstrom model to 7.5% by the integrated ANN model.

Cartesian Genetic Programming evolved Artificial Neural Network (CGPANN)

CGPANN is a Neuro-Evolutionary technique based on Cartesian Genetic Programming where the nodes are replaced by articial neurons, the basic elements of CGPANN. A CGPANN node or neuron is composed of inputs, weighted connections, non-linear activation functions and outputs. The genotype is generated by encoding all these attributes of the neural networks [START_REF] Khan | Efficient representation of recurrent neural networks for markovian/non-markovian non-linear control problems[END_REF]. The parent genotype undergoes mutation, producing the offspring. The network and its parameters continue to evolve until the best possible solution is obtained [START_REF] Miller | Cartesian genetic programming[END_REF]. In this paper, CGPANN is used as the principle estimator in the prediction system. For the production of the next generation population in CGPANN, (1+) evolutionary strategy is recommended. The 10 % mutation rate used implies that only 10% of the genes take part in mutation to give the offspring. During the process of evolution some neurons, whose outputs remain unconnected (Junk Nodes) may not participate at all in generating the system output whereas others (Active Nodes) may play an active role in producing the system output. The weights are generated randomly (ranging from -1 to 1). The output can either be the system inputs or the output of any of the nodes. The neurons in the network are not fully connected. The systems inputs are not connected with all of the neurons in the input layer. The fundamental entities of a typical node in a CGPANN can be observed in Figure 1. The selection of the inputs is from the input array I p , such that

I p = [i p 1, i p 2, i p 3, ....., i p n] (1)
The weighing matrix [array of weights] W g , W g = w g1 , w g2 , w g3 , ....., w gn

for a particular genotype g, consists of random weights ranging from +1 and -1 . The output of a junction in ANN that takes x i as junction input, where the input is multiplied by a randomly alloted weight, w gi can be represented by

O n = N i=1 x i w gi (3)
If the inputs to a node are N, output O nj is

O nj = ξ j (O n ) = ξ j N i=1 x i w gi (4)
here, ξ is the node function. Also, j is defined by If I g is the input set to G k , where each input i has unique value:

i R ∧ 1 ≥ i ≥ 0 (6)
then, network G k consists of randomly chosen inputs I g , node outputs, O nj , for a single output O p such that

O p = 1 n N i=1 O i (7)
where

O i = ξ (O nj W gj + O nj-1 W gj-1 + .... + O n1 W g1 + IW gk ) (8)
Here, W gj , W gj-1 , ... W g1 are Random subsets of W such that W gk is a subset W gj and

W gj = {w gk |w gk R ∧ +1 ≥ w k ≥ -1} (9)
now

G k = {I, O nj , O nj-1 , ...., O n1 , O p } (10)
if we have two successive genotypes G k and G l . G l is the product of G k by mutating µ% characteristics of the genotypes. Let the total entries in G k be

N cf w then N cf w = µ × N cf w (11)
N c f w represents genotype entries that are to be mutated to get G l . The set ζ has ω as a unique entry such that

{ω|ω ζ∀(W g , I, O nj )} (12)
Each value in ω is defined Pseudo-randomly. N c f w values are provided by the available N c f w entities using the equivalent expression ω i |ω i {1, 2, 3, ..., N cf w } ∧ ω i ⊂ {1, 2, 3, ..., N cf w } , i = 1, 2, 3, ...N cf w [START_REF] Sozen | Solar-energy potential in Turkey[END_REF] Where ω i can its value from {N c f w, . . . 3, 2, 1} so each entry in ζ is

γ(i) = ω i |ω i {1, 2, 3, ...N cf w } ∧ {1 ≤ ω i ≤ N cf w } (14)
The network takes inputs from array I p . The network inputs are fed to the primary layer, the output of which is fed into secondary processing nodes, known as intermediate layer nodes. Activation function ξ, that is employed for this particular CGPANN is Log-sigmoid function, given by Eq. 16.

ξ(x) = 1 1 + e -x (15)
Here x is input to the activation function ξ(x). x is defined by the Eq. 17.

O(c, g, j)

= N i=1 w g (c, i).I(g, c, i) (16)
Whereby N defines the quantity of inputs to the given node j for n = 1, 2, 3... N numbered inputs. K are the network inputs in set k = 1, 2, 3...N...K. g represents a particular genotype in a population and c is the node taking part in the operation. w(i, c) is the randomly assigned weight to each input and I(c, g, i) is defined as a randomly selected input to the node j.

M AP E = 1 N N i=1 L F -L A L A × 100% (18)
F itness = 100% -M AP E (19)
Where L F is the predicted value, L A is the actual value and n is the number of days.

Results and Analysis

The proposed CGPANN model has been evaluated using 24 hours input network for the hourly, daily and weekly prediction of solar irradiance. The network takes 24 entries for the prediction of the next hour, next 12 hours, 24 hours, 2 days and 7 days. While the number of nodes on the network is varied between 50 and 500 with a step size of 50, training and testing experiments are repeatedly performed. Subsequently, a phenotype was being translated by the genotype that showed the best fitness level and then testing was carried on hourly, daily, and weekly basis. All the network nodes however may not be utilized by the final phenotype. Usually 5 to 10 percent of nodes participate in the production of the phenotype.

During testing, an independent dataset was given to the network for its validation. Table 1 contains the MAPE values calculated using 24 hours input network for the prediction of the next one hour, 12hrs, and 1 day. The network used here has been trained for hourly prediction. The values change with the changing number of nodes. MAPE values increase as the network is being used for long-term forecasting. The best result has been for the hourly prediction for 500 nodes i.e. a MAPE of only 4.52%.

Table 2 shows the testing results of predictions on 12-hourly, daily, and weekly basis, using the network model that has been trained for 24-hourly prediction. A lowest MAPE value of 5.36 has been achieved for the 12-hourly prediction. Hence, with the increase in time span, the forecasting of the model is affected owing to environmental, technical and climatic issues. The difference in values of Tables 1 and2 explain the fact that no matter what the number of nodes, the CGPANN model exhibits better results when forecasting over a shorter period. However, the model that is trained to predict a single instance is also capable of predicting up to 24 hours without much error (see table 1). And the model that is trained to predict 24 hours can predict up to 7 days (168 hours) in advance with more than 90% accuracy (see table 2), indicating the robustness of the model. Another fact revealed is that the duration for which the network is trained is directly proportional to the efficiency of the model. This is indicated by the lower MAPE values in

Conclusion

The use of renewable energy necessitates exact estimation of renewable potential throughout the globe. Our research is an effort to introduce a technique for solar irradiance forecasting based on the Neuro-evolutionary technique called Cartesian Genetic Programming evolved Artificial Neural Network (CGPANN). The use of this technique could be very useful for improving efficiencies of solar based applications, the integration of renewable energy resources or the management of solar energy systems. The experiments display promising results, yielding an accuracy of 95.48 % and MAPE of 4.52% for hourly prediction. As the only input to the learning model is the solar radiation patterns, in the future, other parameters such as environmental conditions can be incorporated as inputs to produce even more accurate results.

Fig. 3: The actual and the estimated normalized solar irradiance for 300 hours

Fig. 1 :

 1 Fig. 1: The Nth node of CGPANN

IFig. 2 :

 2 Fig. 2: A typical CGPANN genotype and phenotype

Figure 3

 3 Figure3represents a graph that displays the deviation between the actual and the estimated data in terms of the normalized values of irradiance. It is obvious that the estimated data is following the frequent fluctuations in the actual data throughout the time frame i.e. 300 hours. The graph indicates only a minor difference between the two sets of values, confirming the accuracy and robustness of the CGPANN model.

  Ahsa lies at an altitude of 178 m, 25.3 o N latitude and 49.48 o N longitude. The average temperature is 304 K, relative humidity 23% and the daily average solar global radiation is 21.6MJ/m 2 . Although the network has been tested and trained using the data from 1997, the datasets for both belong to different months which mean that the data used for training is not being used for testing. The Mean Absolute Percentage Error (MAPE) value is used as the performance criterion to compare the offspring to the fittest network (the parent) and decide which of these are fit enough to be promoted to the next generation. The parent network further produces nine more networks by mutation. The process continues until the desired fitness is achieved or the maximum number of generations is completed. For training the network, all experiments are run for one million generations. Because a mutation rate of 10% is thought to produce more satisfactory results in comparison to other rates, it has been used in our experiment(Kadilar & Alada, 2009; Huang et al., 2006; Chen et al., 2008).MAPE is given by:

3.1 Experimental Setup

The solar irradiance data of 1997, collected from Saudi Arabian stations in Al-Ahsa

[26]

, were used to train the network for hourly and 24-hourly prediction and subsequently testing was done on the data acquired from 1997 as well as 1998, for hourly, 12-, 24-, 48-hourly and 7 days prediction capability of the model. The inputs to the network are the solar irradiance values, where irradiance, given in W/m2 , is the intensity of electromagnetic radiation incident per unit area. Al-

  Table 2 corresponding to 24-hourly prediction training compared to those in Table 1 that correspond to hourly prediction training.

Table 1 :

 1 Testing Results for Different Number of Nodes using 24 Hours Input Network for Prediction of the Next One Hour, 12 hours, and 1 Day, for Hourly Prediction Training

	Prediction Year 50 100 150 200 250 300 350 400 450 500
	Hourly	1998 6.65 4.83 4.61 4.79 4.65 4.77 5.05 5.15 4.94 4.52
	12-hourly 1997 9.88 7.48 7.03 7.06 6.98 6.64 7.44 7.57 7.22 6.78
		1998 11.01 8.37 7.97 7.96 7.87 7.55 8.36 8.49 8.14 7.68
	24-hourly 1997 10.48 8.42 7.99 8.20 8.09 7.55 8.73 8.94 8.31 7.83
		1998 11.72 9.50 9.10 9.27 9.13 8.62 9.83 10.03 9.41 8.88

Table 2 :

 2 Testing Results for various Number of Nodes for 24 Hours Input Network and Prediction of the Next 12 hours, 1 day, 2 Days, and 1 Week, for 24-hourly Prediction Training Prediction Year 50 100 150 200 50 300 350 400 450 500 12-hourly 1997 6.78 5.41 5.37 5.36 5.39 5.42 5.41 5.38 5.42 5.37 1998 8.15 6.67 6.64 6.63 6.65 6.67 6.67 6.65 6.68 6.64 24-hourly 1998 8.08 6.68 6.65 6.64 6.66 6.70 6.68 6.66 6.69 6.65 48-hourly 1997 7.80 5.98 5.92 5.90 5.95 6.03 5.99 5.93 6.00 5.92 1998 9.41 7.45 7.41 7.41 7.43 7.50 7.45 7.42 7.46 7.42 7 days 1997 8.30 7.43 7.28 7.17 7.39 7.52 7.56 7.26 7.59 7.31 1998 12.41 8.52 8.37 8.29 8.48 8.66 8.62 8.36 8.66 8.40


Engin Yesil
email: yesileng@itu.edu.tr

Tufan Kumbasar
email: kumbasart@itu.edu.tr

M Furkan Dodurka

Ahmet Sakalli
email: sakallia@itu.edu.tr

Peak Observer based Self-tuning of Type-2 Fuzzy PID Controllers

Keywords: Type-2 Fuzzy PID controllers, self-tuning, peak observer

Fuzzy PID (proportional-integral-derivative) controllers are commonly used as an alternative to the conventional PID controllers. In order to improve the control system performance of these controllers many self-tuning methods are already studied. It is mostly observed that the self-tuning mechanism should tune the scaling factors of the fuzzy controller to enhance the transient system performance. On the other hand, these studies only focus on the ordinary (Type-1) Fuzzy PID controllers. In this study, Type-2 Fuzzy PID controllers are studied and a peak observer based self-tuning method is proposed for these controllers. In order to show the benefit of the proposed approach, several Matlab simulations are performed where different type of fuzzy control structures are compared. The results obtained from the simulation studies clearly show the advantage of the proposed approach.

Introduction

Ordinary (type-1) Fuzzy PID (proportional-integral-derivative) controllers (FPID) are accepted as an alternative to conventional PID controllers since they are analogous to the PID controllers from the input-output relationship point of view [START_REF] Galichet | Fuzzy controllers: synthesis and equivalences[END_REF][START_REF] Huang | A fuzzy PID controller being like parameter varying PID[END_REF][START_REF] Qiao | PID type fuzzy controller and parameters adaptive method[END_REF]. Numerous techniques have been developed in the literature for analyzing and designing a wide variety of Type-1 Fuzzy PID (T1-FPID) control systems [START_REF] Qiao | PID type fuzzy controller and parameters adaptive method[END_REF][START_REF] Li | Conventional fuzzy control and its enhancement[END_REF][START_REF] Duan | Effective tuning method for fuzzy PID with internal model control[END_REF][START_REF] Hu | New methodology for analytical and optimal design of fuzzy PID controllers[END_REF]. The design parameters of the T1-FPID controllers can be summarized within two groups, structural parameters and tuning parameters [START_REF] Hu | New methodology for analytical and optimal design of fuzzy PID controllers[END_REF]. The structural parameters include input/output variables to fuzzy inference such as fuzzy sets, Membership Functions (MFs) shapes, rules and inference mechanism. Tuning parameters include input/output Scaling Factors (SFs) and the parameters of the MFs. Starting with the pioneering study by Qiao and Mizumoto [START_REF] Qiao | PID type fuzzy controller and parameters adaptive method[END_REF], various online tuning mechanisms have been presented to improve the control performance of the fuzzy control system in presence of parameter variations and nonlinearities [START_REF] Mudi | A robust self-tuning scheme for PI-and PD-type fuzzy controllers[END_REF][START_REF] Woo | A PID-type fuzzy controller with self-tuning scaling factors[END_REF][START_REF] Guzelkaya | Self-tuning of PID-type fuzzy logic controller coefficients via relative rate observer[END_REF][START_REF] Ahn | Online tuning fuzzy PID controller using robust extended Kalman filter[END_REF]. Thus, Self-Tuning Type-1 Fuzzy PID (STT1-FPID) controllers have been proposed where the SFs or the parameters of MFs have been adjusted in an online manner [START_REF] Mudi | A robust self-tuning scheme for PI-and PD-type fuzzy controllers[END_REF][START_REF] Woo | A PID-type fuzzy controller with self-tuning scaling factors[END_REF][START_REF] Guzelkaya | Self-tuning of PID-type fuzzy logic controller coefficients via relative rate observer[END_REF][START_REF] Ahn | Online tuning fuzzy PID controller using robust extended Kalman filter[END_REF][START_REF] Karasakal | An error-based on-line rule weight adjustment method for fuzzy PID controllers[END_REF][START_REF] Karasakal | Online tuning of fuzzy PID controllers via rule weighting based on normalized acceleration[END_REF]. However, the main research focus was on the tuning of the SFs since their effect on the system's response can easily be observed [START_REF] Woo | A PID-type fuzzy controller with self-tuning scaling factors[END_REF][START_REF] Guzelkaya | Self-tuning of PID-type fuzzy logic controller coefficients via relative rate observer[END_REF]. The presented self-tuning mechanisms for the SFs provide extra degrees of freedom to the fuzzy control structure, and also enhance new tuning structures which have to be determined as updating functions or extra fuzzy inference mechanisms.

Recently, researchers began investigating Interval Type-2 Fuzzy Logic Controllers (IT2-FLCs) which have demonstrated significant control performance improvements in comparison to its type-1 counterpart [START_REF] Wu | Interval Type-2 Fuzzy PI Controllers: Why They Are More Robust[END_REF][START_REF] Wu | On the Fundamental Differences between Type-1 and Interval Type-2 Fuzzy Logic Controllers[END_REF][START_REF] Hagras | A Hierarchical Type-2 Fuzzy Logic Control Architecture for Autonomous Mobile Robots[END_REF]. It has been shown in various works that the T1-FPID controllers using Type-1 Fuzzy Sets (T1-FSs) might not be able to fully handle the high levels of uncertainties associated with control applications while the Interval Type-2 Fuzzy PID (IT2-FPID) controller using Interval Type-2 Fuzzy Sets (IT2-FSs) might be able to handle such uncertainties to produce a better control performance [START_REF] Hagras | A Hierarchical Type-2 Fuzzy Logic Control Architecture for Autonomous Mobile Robots[END_REF][START_REF] Castillo | Type-2 Fuzzy Logic Theory and Applications[END_REF][START_REF] Wu | A simplified type-2 fuzzy logic controller for real-time control[END_REF]. It has been shown that IT2-FPIDs achieve better control performances because of the additional degree of freedom provided by the Footprint of Uncertainty (FOU) in their antecedent MFs [START_REF] Wu | On the Fundamental Differences between Type-1 and Interval Type-2 Fuzzy Logic Controllers[END_REF][START_REF] Hagras | A Hierarchical Type-2 Fuzzy Logic Control Architecture for Autonomous Mobile Robots[END_REF][START_REF] Kumbasar | Type-2 fuzzy model based controller design for neutralization processes[END_REF]. In literature, the design of the IT2-FPID controllers is usually solved by extending MFs of an existing T1-FPID controller or by employing optimization algorithms [START_REF] Kumbasar | Big Bang-Big Crunch optimization based interval type-2 fuzzy PID cascade controller design strategy[END_REF][START_REF] Kumbasar | A simple design method for interval type-2 fuzzy PID controllers[END_REF][START_REF] Yesil | Interval type-2 fuzzy PID load frequency controller using Big Bang-Big Crunch optimization[END_REF].

In this study, the peak observer (PO) based self-tuning method, which was first proposed for T1-FPID controllers, is extended for IT2-FPID controllers. Till now, there is a very limited research on the self-tuning methods for the Type-2 Fuzzy Logic Controllers. The benefit of this approach is examined via various simulation studies performed on First Order plus Dead Time (FOPDT) systems. The simulation results clearly show that IT2-FPID controller achieve a better performance than T1-FPID and PO based T1-FPID controllers. When the proposed self-tuning approach is employed to the IT2-FPID controller the control system performance increases significantly.

The Design of Fuzzy PID Controllers

Interval Type-2 Fuzzy PID controllers are generalized forms of Type-1 Fuzzy PID controllers, which are mostly mentioned as an alternative to conventional PID controllers since they are analogous to the PID controllers from the input-output relationship point of view [START_REF] Galichet | Fuzzy controllers: synthesis and equivalences[END_REF][START_REF] Huang | A fuzzy PID controller being like parameter varying PID[END_REF][START_REF] Qiao | PID type fuzzy controller and parameters adaptive method[END_REF]. PID controllers are widely used in many control areas such as process control, adaptive control, robust control, nonlinear control, automation and robotics and they are still the most used controller type in the industry [START_REF] Astrom | Advanced PID Control[END_REF]. The design strategies of conventional PID controllers depends on the mathematical model of the system, however conventional PID controllers may not provide satisfactory results in case of imperfect model. Conventional PID controllers provide a linear transformation or mapping between their inputs and outputs while Fuzzy PID controllers provide a nonlinear transformation or mapping between their inputs and outputs.

In literature, various fuzzy controller structures are proposed including PD-type (Proportional-Derivative), PI-type (Proportional-Integral) and PID-type [START_REF] Li | Conventional fuzzy control and its enhancement[END_REF][START_REF] Guzelkaya | Self-tuning of PID-type fuzzy logic controller coefficients via relative rate observer[END_REF]. In this context, fuzzy controllers are constructed by choosing the inputs to the error (e) and derivative of the error (𝑒) and the output as the control signal (u). The inputs of Fuzzy Logic Controllers (FLCs) (regardless to the type of the controller) is as follow

𝑒 = 𝑟 -𝑦 (1)
𝑒̇= 𝑑𝑒 𝑑𝑡 = 𝑑(𝑟 -𝑦) 𝑑𝑡 ( 2
)
where 𝑟 is the reference signal or setpoint and 𝑦 is the system output. In fuzzy controller structures (regardless to the type of the controller), input SFs normalize the inputs to the universe of discourse, particularly, 𝐾 𝑒 normalizes the error signal and 𝐾 𝑑 normalize the derivative of error signal as follow

𝐸 = 𝐾 𝑒 𝑒, 𝐸 ∈ [-1,1] (2)
𝐸 ̇= 𝐾 𝑑 𝑒, 𝐸 ̇∈ [-1,1] (3)
where 𝐸 is the normalized error signal of the FLC and 𝐸 ̇is the normalized derivative of error signal of the FLC. Here, the universe of discourses of the inputs (𝐸 and 𝐸 ̇) is defined in the interval [-1, 1]. The FPID controller structure including inputs, SFs and the FLC (T1-FLC or IT2-FLC) is illustrated as in Fig. 1. 𝐾 𝑎 is a constant gain, and 𝐾 𝑏 is the integral gain of the output of the FLC. Here also notice that, output SFs make the control signal applicable for a particular problem or a real life problem because the output of the FLC could only take values in the interval [-1,1] due to its universe of discourse of the output MFs. Then the generalized formulation of the output of a FLC that used in this study is as follows:

𝑢 = 𝐾 𝑎 𝑈 + 𝐾 𝑏 ∫ 𝑈 𝑑𝑡, 𝑈 ∈ [-1,1] (4)

The Type-1 Fuzzy PID Controller Structure

In this subsection, the internal structure and the design strategies of the employed T1-FPID controllers are presented. The rule structure is as follows:

𝑅 𝑞 : If 𝐸 is 𝐴 1𝑗 and 𝐸 ̇ is 𝐴 2𝑗 THEN 𝑈 is 𝐶 𝑞 , 𝑖, 𝑗 = 1,2,3; 𝑞 = 1, … , 𝑄 = 9 (5)
where 𝐴 1𝑗 is and 𝐴 2𝑗 are the antecedent MFs for the inputs 𝐸 and 𝐸 ̇ , respectively, 𝐶 𝑞 is the consequent crisp set and 𝑄 is the number of rules. In the handled T1-FPID structure, a symmetrical 3x3 rule base is used as shown in Fig. 2a. Here, three Type-1 Fuzzy Sets (T1-FSs) for each input domain (𝐸 and 𝐸 ̇) are used and they are denoted as N (Negative), Z (Zero) and P (Positive). The T1-FSs of the T1-FLC are defined in [-1 1] interval as shown in Fig. 2b. The consequent part is defined with five singleton consequents, as shown in Fig. 2c. The implemented T1-FLCs use the product implication and the weighted average defuzzification method.

The Interval Type-2 Fuzzy PID Controller Structure

In this subsection, the general structure of the employed IT2-FPID controller is presented. In the design of the IT2-FPID, the same symmetrical 3x3 rule base given in Fig. 2a and consequent MFs of T1-FPID illustrated in Fig. 2c are used. Thus, the rule structure of the IT2-FLC is as follows:

𝑅 𝑞 : IF 𝐸 𝑖𝑠 𝐴 ̃1𝑖 and 𝐸 ̇ 𝑖𝑠 𝐴 ̃2𝑗 THEN 𝑈 is 𝐶 𝑞 , 𝑖, 𝑗 = 1,2,3; 𝑞 = 1, … , 𝑄 = 9 (6)
where 𝐶 𝑞 is the consequent crisp set, 𝑄 is the number of rules and 𝐴 ̃1𝑖 and 𝐴 ̃2𝑗 are the antecedent Interval Type-2 Fuzzy Sets (IT2-FSs). We will denote the IT2-FSs as N (Negative), Z (Zero) and P (Positive). The antecedent IT2-FSs can be described in terms of upper MFs (𝜇 𝐴 ̃1𝑖 and 𝜇 𝐴 ̃2𝑖

) and lower MFs (𝜇 𝐴 ̃1𝑖 and 𝜇 𝐴 ̃2𝑖 ) to form the Foot of Uncertainty (FOU), which provides extra degree of freedom in IT2-FSs. The antecedent IT2-FSs are defined again in [-1 1] interval. Thus, the FOU will be created by the heights of the lower MFs (m1 and m2) of the IT2-FLC. The implemented IT2-FLC uses the center of sets type reduction/ defuzzification method [START_REF] Liang | Interval type-2 fuzzy logic systems: theory and design[END_REF]. The Peak Observer based Self-tuning Method

As shown in Fig. 1, the actual control signal of the closed loop system is a function of scaling factors and the output of the FLC as 𝑢 = 𝑓(𝑈, 𝐾 𝑎 , 𝐾 𝑏 ) or 𝑢 = 𝑓(𝑈, ∫ 𝑈 , 𝐾 𝑎 , 𝐾 𝑏 ). The FLC performs a nonlinear mapping between its inputs and output, so the output of the FLC can also be expressed as a function of its inputs namely, error and derivative of the error as 𝑈 = 𝑓(𝐸, 𝐸 ̇) or 𝑈 = 𝑓(𝑒, 𝑒, 𝐾 𝑒 , 𝐾 𝑑 ). Note that 𝑓(. ) represents a functional mapping. Consequently, the FPID controller structure becomes analogous to conventional PID structure. This analysis was examined in detail by considering the input space into several subspaces and mathematical expression of MFs [START_REF] Qiao | PID type fuzzy controller and parameters adaptive method[END_REF] and the control signal can be written as 𝑢 = 𝐾 𝑎 𝐾 0 + 𝐾 𝑏 𝐾 0 𝑡 + (𝐾 𝑎 𝐾 𝑒 𝐾 1 + 𝐾 𝑏 𝐾 𝑑 𝐾 2 )𝑒 + (𝐾 𝑏 𝐾 𝑒 𝐾 1 ) ∫ 𝑒 𝑑𝑡 + (𝐾 𝑎 𝐾 𝑑 𝐾 2 )𝑒̇ [START_REF] Mudi | A robust self-tuning scheme for PI-and PD-type fuzzy controllers[END_REF] where 𝐾 0 , 𝐾 1 and 𝐾 2 are the constants which are calculated with respect to antecendent and consequent membership functions. They are calculated with respect to values of consequent sets and apexes of the antecedent membership functions for corresponding inputs and subspaces. Here, 𝐾 0 is a nonlinear term, which converges to zero by time. So one can say that a FPID controller is similar to a conventional PID controller with variable gains. Then after neglecting the nonlinear term for the steady state, the PID equivalent control components of FPID can be achieved as follows

𝐾 𝑝 = 𝐾 𝑎 𝐾 𝑒 𝐾 1 + 𝐾 𝑏 𝐾 𝑑 𝐾 2 (8)
𝐾 𝐼 = 𝐾 𝑏 𝐾 𝑒 𝐾 1

𝐾 𝐷 = 𝐾 𝑎 𝐾 𝑑 𝐾 2 [START_REF] Ahn | Online tuning fuzzy PID controller using robust extended Kalman filter[END_REF] where Kp is the proportional gain, KI is the integral gain and the KD is the derivative gain of the fuzzy controller.

Fig. 4. Peak observer based self-tuning for Fuzzy PID controllers

The peak observer based self-tuning method for T1-FPID controllers has been proposed in [START_REF] Qiao | PID type fuzzy controller and parameters adaptive method[END_REF]. The block diagram of the proposed method is shown in Fig. 4. The peak observer observes the output of the system and transmits a signal at each peak time and measures the absolute peak value. The parameter regulator tunes the controller parameters Kd and Kb simultaneously at each peak time according to the peak value. The algorithm for tuning these parameters is as follows:

𝐾 𝑑 =
𝐾 𝑑 𝑠 𝛿 𝑘 , 𝐾 𝑏 = 𝛿 𝑘 𝐾 𝑏 𝑠 [START_REF] Karasakal | An error-based on-line rule weight adjustment method for fuzzy PID controllers[END_REF] where 𝐾 𝑑 𝑠 and 𝐾 𝑏 𝑠 are the initial values of 𝐾 𝑑 and 𝐾 𝑏 , 𝛿 𝑘 are the peak values. It can easily be deduced from the relation given in [START_REF] Guzelkaya | Self-tuning of PID-type fuzzy logic controller coefficients via relative rate observer[END_REF] that if in the meanwhile of decreasing 𝐾 𝑏 , 𝐾 𝑑 is increased in the same rate as 𝐾 𝑏 is decreased, the equivalent proportional control strength will remain unchanged. Then, the system can always keep quick reaction against the error under this condition. This is achieved by updating the integral coefficient as the reciprocal of the derivative coefficient. As it is seen from ( 12), the integral factor is decreased in a direct proportion to the peak error value.

Simulation Studies

In this section, the performance of the proposed self-tuning method for IT2-FPID controllers is examined via various cooperative simulations. For this purpose, different processes are used for simulations but as it is well-known a large number of industrial plants can be modeled by a first-order plus dead time (FOPDT) transfer function. The wide use of FOPDT is due both to its simplicity as well as its ability to capture the essential dynamics of several industrial processes. Therefore, the results obtained only for FOPDT process are presented here. The transfer function of FOPDT is as follows:

𝐺(𝑠) = 𝐾 𝑇𝑠 + 1 𝑒 -𝐿𝑠 ( 12
)
where K is the process gain, T is the time constant and finally L is the time delay.

In the simulations, four different control schemes, namely Type-1 Fuzzy PID Controller (T1-FPID), Peak Observer based Type-1 Fuzzy PID Controller (POT1-FPID), Interval Type-2 Fuzzy PID Controller (IT2-FPID), and the proposed Peak Observer based Interval Type-2 Fuzzy PID Controller (POIT2-FPID), are used. In order to compare the performance of transient responses of the proposed POIT2-FPID controller with other fuzzy controllers, the following four performance measures are considered: i) Integral Square Error (ISE) which is defined as

  dt ) t ( y ) t ( r ISE 0 2     (13)
ii) Integral Absolute Error (IAE) which is defined as

dt ) t ( y ) t ( r IAE 0     (14)
iii) Integral Time Absolute Error (ITAE) which is defined as

dt ) t ( y ) t ( r t ITAE 0     (15)
iv) Total Variation (TV) [START_REF] Karasakal | An error-based on-line rule weight adjustment method for fuzzy PID controllers[END_REF] of the control input u(t), which is defined as

      1 i i 1 i u u TV (16)
As the first step, the nominal process parameters are set to K=1, T=1 and L=0.2. The sampling time of the simulation is set to 0.05 s and noise is added to the process output to make the simulations more realistic. Then the input SFs of the T1-FPID controller are designed in order to have a response without an overshoot but small rise time and settling time. For this purpose, first the input scaling factor of the error (Ke) is calculated as

𝐾 𝑒 = 1 𝑟(𝑡 𝑓 ) -𝑦(𝑡 𝑓 ) (17)
where 𝑟(𝑡 𝑓 )and 𝑦(𝑡 𝑓 ) are the values of the reference and system output at the time of the reference variation. Since the unit step changes are studied for the simulations Ke is set to 1, then the other scaling factors are determined as Kd = 0.5141, Ka =0.077, Kb = 7.336 after optimization. As the second step, the IT2-FPID controller is design using the same scaling factors obtained for T1-FPID controller. For the design procedure only the height of the lower membership functions (m1 and m2) are designed. After the optimization procedure, m1 and m2 are set to 0.3 and 0.9, respectively. The step responses of the process for four different control schemes are illustrated in Fig. 5. As clearly seen from Fig. 5, the set-point following performance of all fuzzy controllers are very close to each other for the nominal process. Since there is no overshoot, the self-tuning mechanism did not tune the SFs. In order to show the benefit of using IT2-PID controllers, and the proposed PO based Type-2 Fuzzy PID controller the system parameters of the process are perturbed. The step responses of four controller for the first perturbed process, which has the parameters as K=1.3, T=1.9, L=0.4, is given in Fig. 6a. For this case, IT2-FPID controller gives a better set-point following performance than the T1-FPID controller because of its extra degree of freedom. As it is seen in Fig. 6b and Fig. 6c, the two scaling factors (Kd and Kb) of the fuzzy controllers are tuned as the overshoots occurs, therefore the performance of the fuzzy controllers is improved drastically. The control signals illustrated in Fig. 6d, shows that IT2-FPID controller with peak-observer based self-tuning mechanism has a smoother control action in comparison to the other fuzzy controllers. The simulation results obtained for the first perturbed process show that IT2-FPID controller has a better set-point following performance than the T1-FPID controller because of FOU. Moreover, the peak observer based self-tuning mechanism improves the performance of the fuzzy controllers. The best performance from the point of transient response is obtained by the proposed IT2-FPID controller with the peak observer mechanism. Moreover, two more simulations are performed for two other perturbed processes and the control performances are given in Fig. 7. Parallel to the first simulation results, IT2-FPID controller always outperforms the T1-FPID controller. The Peak Observer based T1-FPID is more robust against parametric uncertainties; consequently the process output inherits smaller oscillations. When the proposed Peak Observer selftuning method is employed to IT2-FPID controller, its performance again increases drastically.

The performance results of the fuzzy controllers are tabulated in Table 1. As it is seen from Table 1, the IT2-FPID controller is always better than the T1-FPID controller because of the additional degree of freedom provided by the FOU in their antecedent MFs. Also, the peak observer self-tuning methods improve the fuzzy control system performance. On the other hand, since the peak observer changes the scaling factors of the fuzzy controller, the control signal varies and thus results with higher TV values compared to the counterparts without a self-tuning mechanism.

Conclusion

In this study, the peak observer based self-tuning method, which involves a peak observer and parameter regulator mechanisms, is extended to IT2-FPID Controllers. The advantage of the proposed approach is tested via Matlab simulations and the control system performance of the proposed method is compared with T1-PID controllers with and without a self-tuning mechanism, and also with an IT2-FPID controller without a self-tuning method. The obtained result clearly shows that peak-observer based selftuning approach improved the performance in comparison to the IT2-FPID controller. As a future work, this approach will be implemented in a PLC and real-time experiments will be performed.

Fig. 1 .

 1 Fig. 1. Fuzzy PID Controller

Fig. 2 .

 2 Fig. 2. (a) The fuzzy rule base, (b) the antecedent MFs, (c) the consequent MFs used for the FPID controllers.

Fig. 3 .

 3 Fig. 3. Interval Type-2 Fuzzy Sets used for IT2-PID controller

Fig. 5 .

 5 Fig. 5. The set-point following performance of the fuzzy controllers for the nominal process.

Fig. 6 .

 6 Fig. 6. The first the perturbed process (K=1.3, T=1.9, L=0.4): (a) The process output; (b) The change of Kb; (c) The change of Kd; (d) The control signals.

Fig. 7 .

 7 Fig. 7. The set-point following performance of the fuzzy controllers for the perturbed processes: (a) K=1.1, T=1.3, L=0.45; (b) K=3, T=3, L=0.3.

Table 1 .

 1 Simulation comparison using different performance measures.

	Process	Controller	ISE	IAE	ITAE	TV
		T1-FPID	0.552	0.842	1.564	0.409
	Nominal Process	POT1-FPID	0.547	0.872	1.745	0.970
	K=1, T=1, L=0.2	IT2-FPID	0.550	0.842	1.590	0.364
		POIT2-FPID	0.550	0.838	1.508	0.926
	Perturbed Process K=1.3, T=1.9, L=0.4	T1-FPID POT1-FPID IT2-FPID POIT2-FPID	1.272 1.050 1.160 1.055	2.677 1.755 2.181 1.733	9.430 3.357 5.497 3.240	1.589 1.401 1.000 1.275
	Perturbed Process K=1.1, T=1.3, L=0.45	T1-FPID POT1-FPID IT2-FPID POIT2-FPID	2.969 1.457 1.503 1.252	5.739 2.423 2.981 2.062	36.223 5.659 10.394 4.310	3.756 1.690 1.476 1.461
		T1-FPID	0.836	1.710	4.355	1.166
	Perturbed Process	POT1-FPID	0.789	1.488	3.278	1.290
	K=3, T=3, L=0.3	IT2-FPID	0.789	1.472	2.979	0.768
		POIT2-FPID	0.763	1.382	2.940	1.173

Acknowledgments. This research is supported by the Scientific and Technological Research Council of Turkey (TUBITAK) under the project (113E206). All of these supports are appreciated.