Research Article | Open Access
An Improved Global Harmony Search Algorithm for the Identification of Nonlinear Discrete-Time Systems Based on Volterra Filter Modeling
This paper describes an improved global harmony search (IGHS) algorithm for identifying the nonlinear discrete-time systems based on second-order Volterra model. The IGHS is an improved version of the novel global harmony search (NGHS) algorithm, and it makes two significant improvements on the NGHS. First, the genetic mutation operation is modified by combining normal distribution and Cauchy distribution, which enables the IGHS to fully explore and exploit the solution space. Second, an opposition-based learning (OBL) is introduced and modified to improve the quality of harmony vectors. The IGHS algorithm is implemented on two numerical examples, and they are nonlinear discrete-time rational system and the real heat exchanger, respectively. The results of the IGHS are compared with those of the other three methods, and it has been verified to be more effective than the other three methods on solving the above two problems with different input signals and system memory sizes.
The Volterra model is a kind of nonlinear filter model, which is usually employed to track and identify plenty of complex nonlinear systems. In order to enhance the quality of the system identification, it is very crucial issue to select optimum model coefficient called the kernel. Therefore, the Volterra model is essentially an extension of linear filter model to nonlinear case. During the past decade, there have been many research works on the Volterra model. Campello et al.  tackled the problem of expanding Volterra models using Laguerre functions. The global optimal solution is obtained when each multidimensional kernel of the model is decomposed into a number of independent orthonormal bases. Furthermore, the solution obtained is able to minimize the upper bound of the squared norm of the error resulting from the practical truncation of the Laguerre series expansion into a finite number of functions. Masugi and Takuma  described a Volterra system-based nonlinear study of video-packet transmission over IP networks. Based on the Volterra system, the authors performed a time-series analysis of measured data for network response evaluation. The novel method can reproduce the time-series responses observed in video-packet transmission over the Internet, characterizing nonlinear dynamic behaviors such that the obtained results gave an appropriate depiction of network conditions at different times. Gruber et al.  presented a nonlinear model predictive control (NMPC) method based on a second-order Volterra series model for greenhouse temperature control using natural ventilation. These models, denoting the simple and logical extension of convolution models, are capable of describing the nonlinear dynamic feature of the ventilation and other environmental conditions on the greenhouse temperature. Many applications of Volterra series modeling were executed in the frequency-domain based on the Generalized Frequency Response Functions (GFRF) [4, 5]. In the light of the wide applications of Volterra series, Li and Billings  presented an approach to estimate the GFRF in a piecewise manner for duffing type oscillators for the underrepresented weakly nonlinear region. Additionally, they devised a new method to obtain the energy contributions of each order kernel. The new nonparametric method can not only construct the amplitude-invariant GFRF over a certain excitation range, but also avoid building large sets of time-domain models. In addition, the Volterra filter model can be found in some other application areas [7–14].
Chang  devised an improved particle swarm optimization (IPSO) to implement system identification based on Volterra filter model. In this paper, we develop an improved global harmony search (IGHS) algorithm and try the IGHS as an efficient candidate for system identification based on Volterra filter model. The harmony search (HS) algorithm was firstly proposed by . The HS is a simple but efficient algorithm, and its many improved versions have been applied into many problems including reliability problems , reactor core fuel management optimization , and sizing optimization of truss structures .
The paper is organized as follows. In Section 2, the pruned second-order Volterra filter model is simply presented. In Section 3, a novel global harmony search algorithm is introduced. In Section 4, an improved global harmony search algorithm is proposed, and its procedure is fully explained. In Section 5, four harmony search algorithms are used for two examples with different signal inputs and memory sizes. We end this paper with some conclusions in Section 6.
2. Volterra Filter Model and Its Pruned Second-Order Form
The Volterra filter model is an efficient method for the identification of nonlinear discrete systems, and it has come into researchers’ notice in recent decades. The discrete form of Volterra filter model of the th order  is given bywhere denotes the system memory size. Equation (1) denotes the Volterra filter model with the infinite series. However, this model is hard to compute and master due to its complex and expatiatory formula. In this paper, we only study its simplified and approximate form called the truncated second-order Volterra model [7, 19], which is stated as follows:
In order to facilitate expression, (2) can also be expressed as the following vector form:
Here, the superscript represents the transpose of a vector and stands for the Volterra kernel vector given by
In addition, denotes the Volterra input vector given by
To achieve the nearest approximation of the actual system output, appropriate kernel vector should be determined under the input vector . In this paper, an improved global harmony search (IGHS) algorithm is proposed to determine kernel vector . The IGHS is an improved version of novel global harmony search (NGHS) algorithm ; thus, both the NGHS and the IGHS will be presented in the following sections.
3. Novel Global Harmony Search (NGHS) Algorithm
Novel global harmony search (NGHS) algorithm  is a variant of harmony search (HS) algorithm , and it is superior to the HS for solving unconstrained optimization problems. The NGHS improvises new harmony vectors by combining position updating and mutation. Concretely, the steps of NGHS are explained as follows.
Step 1 (initialize the NGHS parameters and the problem parameters). The NGHS parameters consist of harmony memory size , the number of improvisations , and mutation rate . In addition, the problem parameters include the number of problem variables , the lower bound , and upper bound of the th () problem variable. Furthermore, the number of improvisations () is actually the total number of generations for adjusting the parameters related to the identification problem. In each generation, only one new candidate solution including parameters is generated, and this solution is accepted if and only if it is better than the worst one of the previous solutions.
Step 2 (initialize harmony memory (HM)). The initial harmony memory () can be expressed in the following matrix form:where stands for the th () variable of the th () harmony vector. Moreover, it is randomly produced from a uniform distribution in the ranges . In HM, any vector () represents a candidate solution of the parameters (as in (4)) needed for solving the identification problem. Furthermore, the length of is equal to , which is exactly the number of the parameters in the Volterra kernel vector (as in (6)).
Step 3 (improvise a new harmony). Improvisation is actually the operation of producing a new harmony vector. For the NGHS, its improvisation mainly includes two steps, and they are position updating and genetic mutation, respectively. More specifically, the improvisation can be presented in Table 1.
Here, “best" and “worst” stand for the indexes of the best harmony and the worst harmony in , respectively. , , and denote three uniformly generated random numbers in . With respect to position updating, a new harmony vector is improvised near the best harmony vector, which can facilitate the convergence rate of the NGHS. On the other hand, it is worth noticing that genetic mutation is an event of small probability, and it is utilized to avoid the premature convergence of the NGHS.
Step 4 (update harmony memory). Replace the worst harmony vector of with the new improvised harmony vector no matter whether is better or worse than .
4. An Improved Global Harmony Search (IGHS) Algorithm
In this paper, an improved global harmony search (IGHS) algorithm is proposed to implement the system identification based on Volterra filter model. The IGHS is an improved version of the NGHS, and it is different from the NGHS in the following two aspects.
(1) The Modification of Genetic Mutation. The genetic mutation of the NGHS is conducted according to uniform distribution. To further explore and exploit the solution space, the uniform distribution is replaced with the other two probability distributions, and they are normal distribution [21–23] and Cauchy distribution [22, 23], respectively. Therefore, the new genetic mutation are implemented by using the following two equations:
Here, denotes the index of the th component of the improvised harmony vector, and denotes the th component of the best harmony vector in . stands for the normal distribution with mean and standard deviation , and stands for the Cauchy distribution with location parameter and scale parameter . If overflow happens to generated using either normal distribution or Cauchy distribution, it will be truncated to . On the other hand, if the condition (as in Table 1) is satisfied, either normal distribution or Cauchy distribution will be utilized to carry out genetic mutation, and the probability of using any one of the two distributions is equal to . By adopting normal distribution and Cauchy distribution, the capacity of escaping from the local optimums is enhanced for the IGHS. In the meantime, the solution space can be fully explored and exploited, which is beneficial for improving the quality of the improvised harmony vector .
In addition to the equation of genetic mutation, we also modify the mutation rate () of the NGHS. The parameter is used to determine whether or not a variable of adopts genetic mutation. A large value is beneficial for searching in a large scope, while a small value is helpful to search in a small scope. To well balance the global search and local search of the IGHS, the value decreases dynamically with increasing generations () as follows: where and represent the minimal and maximal mutation rates. () denotes the current generation number, and the IGHS is stopped when the current generation number reaches the maximal generation number . In addition, denotes a fixed generation number, and it is set to here.
(2) Introduce and Modify an Opposition-Based Learning (OBL) . In order to improve the convergence of the IGHS, a method called opposition-based learning (OBL)  is firstly introduced. In short, this technique can be stated as follows:
Here, is the opposite solution of , and is the th variable of . Furthermore, and represent the minimal and the maximal values in the set , and they are expressed as follows:
According to (10), only the opposite value of is employed to produce the corresponding variable in the range of . In other words, once the value is determined, it has an exclusive opposite value which is equal to .
In this paper, we enhance the IGHS improvisation by modifying the OBL technique. More specifically, the value is randomly produced in the range of , and it is given bywhere denotes a random number uniformly generated in . By using (12), the IGHS can yield numerous possible values in the range of . Compared to the OBL technique, its modified version has wider searching space, which is in favor of fine-tuning of the best harmony vectors.
Here, stands for the th () harmony vector in , and denotes its th () variable. After performing the modified OBL operation, each is compared with the corresponding harmony vector . More precisely, if is better than , should be replaced with . It is worth noticing that the modified OBL is an event of small probability. In other words, this technique is an auxiliary step of the IGHS improvisation, and it is mainly used to improve the quality of the improvised harmony vector.
Based on the above detailed illustrations about the two modifications of the NGHS improvisation, the IGHS improvisation can be summarized in Table 2.
By utilizing the IGHS, the identification diagram of nonlinear discrete-time systems based on the pruned second-order Volterra model is shown in Figure 1. Besides, the nomenclatures appearing in this diagram are explained in Table 3.
The goal of the IGHS is to find the optimal kernel vector for Volterra filter model so that the difference between the estimated output and the actual system output is minimized. Thus, it is advisable to find an objective function to satisfy the design requirement. Here we adopt the following objective function :Here, stands for the mean square error (MSE), and denotes the sampling number. By using the IGHS to optimize MSE, we can obtain the most appropriate kernel vector for the second-order Volterra model.
5. Experimental Results and Analysis
To verify the validity of the IGHS on identifying nonlinear system based on second-order Volterra filter model, two examples including the highly nonlinear discrete-time rational system and the real heat exchanger are considered. More specifically, these two examples are explained as follows.
5.1. Example 1
The first example is the highly nonlinear discrete-time rational system , and its mathematical model is stated as follows:
Here, two types of input signals  are considered. The first is defined as Example a whose input is a random signal uniformly produced in . The second is defined as Example b whose input is a cosine signal . For both Examples a and b, the measurement noise is always supposed to be a Gaussian noise of .
In this experiment, the IGHS is used for Example a with . Moreover, the IGHS parameters are set as follows: the maximal mutation rate , the minimal mutation rate , harmony memory size , the number of improvisations , and sampling number . Moreover, is actually the total number of output samples (actual system output or estimated output). According to (14), there are sampling values for actual system output (). After implementing the IGHS, we can obtain the comparison between the estimated output and the actual output in Figure 2.
From Figure 2, it can be seen that the estimated output approximates the actual output well. In addition, the minimal mean square error (MSE) obtained by the IGHS is equal to , which is a satisfactory result. Besides, we use Example a with to investigate the performance of the IGHS for second-order Volterra filter model with large memory size. The parameters of the IGHS are the same as those used for Example a () except , and it is set to in this experiment. By performing the IGHS, we can obtain the comparison between the estimated output and the actual output in Figure 3.
It is clear from Figure 3 that a satisfactory approximation can be attained by utilizing the IGHS. Additionally, the minimal MSE yielded by the IGHS is equal to for Example a with , which provides better modeling capacity.
In addition to random signal, another testing input signal is used to investigate the performance of second-order Volterra filter model using the IGHS. For Example b with and , the IGHS parameters are the same as those for Example a, and Figures 4 and 5 display the comparisons of results for Example b with and . As expected, the difference between the estimated output and the actual output is small in each case. Additionally, two satisfying MSEs can be obtained, respectively, for and by carrying out the IGHS, and they are equal to and , respectively.
5.2. Example 2
The second example is the real heat exchanger and its mathematical model [7, 25] is given by where is the process input denoting the flow rate and is constrained by the range of , stands for the static nonlinearity, and stands for the process output temperature. By combining (16) and (17), the simplified mathematical model can be expressed as the following difference equation:
Here, two types of input signals [7, 25] are considered. The first is defined as Example a whose input is randomly generated in the range . The second is defined as Example b whose input is a sine signal . According to [7, 25], the measurement noise is excluded from the problem model; thus, we have .
In this experiment, the memory size is set to for second-order Volterra filter model. To fulfill the physical input requirement, the modeling input is randomly generated in the range and the testing input is then set to . Moreover, the IGHS parameters used for Example are the same as those of Example . By implementing the IGHS, the outputs of second-order Volterra model are shown in Figures 6 and 7, respectively, for the modeling input and testing input. Based on the careful observations on experimental results, we can confirm that the satisfactory approximation is attained in each case.
5.3. Comparison of the IGHS with the Other Three HSs
In this paper, four methods consisting of the HS , IHS , NGHS , and IGHS are used to solve the above problems with different input signals and system memory sizes. The parameters of these four methods are set as follows. For the HS, harmony memory considering rate , pitch adjusting rate , and bandwidth . For the IHS, , the minimal pitch adjusting rate , the maximal pitch adjusting rate , the minimal bandwidth , and the maximal bandwidth . For the NGHS, mutation rate . For the IGHS, the minimal mutation rate , the maximal mutation rate , and . Note that all the methods are executed under identical conditions (harmony memory size and the number of improvisations). More specifically, the harmony memory size () of each method is set to , and the (Number of improvisations) value of each method is set to and , respectively, for the system memory sizes and . With respect to second-order Volterra filter model, sampling number is set to , and the range of each variable of the kernel vector is set as . Matlab 7.0 is used to execute the above procedure under the environment of Intel(R) Core(TM) i5-2410M CPU @ 2.30 GHz. 20 independent runs are performed in each case, and the optimization results are presented in Table 4.