# WH-EA: An Evolutionary Algorithm for Wiener-Hammerstein System Identification

Authors: J. Zambrano; J. Sanchis; J. M. Herrero; M. Martínez
Journal: Complexity
 (2018)

---

## Abstract

Current methods to identify Wiener-Hammerstein systems using Best Linear Approximation (BLA) involve at least two steps. First, BLA is divided into obtaining front and back linear dynamics of the Wiener-Hammerstein model. Second, a refitting procedure of all parameters is carried out to reduce modelling errors. In this paper, a novel approach to identify Wiener-Hammerstein systems in a single step is proposed. This approach is based on a customized evolutionary algorithm (WH-EA) able to look for the best BLA split, capturing at the same time the process static nonlinearity with high precision. Furthermore, to correct possible errors in BLA estimation, the locations of poles and zeros are subtly modified within an adequate search space to allow a fine-tuning of the model. The performance of the proposed approach is analysed by using a demonstration example and a nonlinear system identification benchmark.


---

## Body

## 1. Introduction

Nonlinearities are present to a greater or lesser extent in all real processes. When nonlinearities are weak, linear models can be successfully used to forecast the evolution of variables or to design control schemes. Currently, a lot of methods to build linear models can be found in the literature [1–5]. However, when nonlinearities are hard, linear models just can be used in a specific operation range. If process operating range is large, cause-effect relationship should be represented by a nonlinear model. An alternative to nonlinear system modelling is to describe the process phenomena using rigorous first-principles formulation [6–8]; nevertheless, in most cases, it can be a very challenging task. Another alternative is the use of soft computing methods for process identification. In this framework, nonlinear system identification has attracted considerable interest of researchers over the past few years. Nowadays, nonlinear identification is an open research topic where some benchmark problems have been proposed [9–12] and real measurement data are available for testing and validate different nonlinear identification methods (e.g., DaISy database [13]).One of the most challenging problems regarding nonlinear system identification is the selection of a good model structure. Currently there are several structures based on neural networks [14], block-oriented models [15, 16], Volterra series [17], NARMAX models [18], and fuzzy models [19], among others. A review of black box methods to nonlinear identification can be found in Suykens and Vandewalle [20].In this paper, block-oriented models are considered, which are a class of nonlinear representations consisting of linear time-invariant (LTI) systems coupled with nonlinear static functions (NL) [21]. Within this class of models, the most popular ones are Wiener (LTI-NL), Hammerstein (NL-LTI), Wiener-Hammerstein (LTI-NL-LTI), and Hammerstein-Wiener (NL-LTI-NL) models [15]. Nowadays, several methods to identify these models can be found in the literature. An interesting classification of contributions that have been developed until the last decade can be found in Lopes dos Santos et al. [22].Block-oriented models are attractive for their simplicity and great capabilities to model nonlinear dynamic systems [23–28]. Specifically, Wiener-Hammerstein models have proved to be able describe several systems like a paralyzed skeletal muscle [29, 30], a limb reflex control system [31], a DC-DC converter [32], a heat exchanger system and a superheater-desuperheater in a boiler system [33], and a thermal process [34], among others [35]. Not only does the study of block-oriented models address parameters estimation, but also these structures are used to implement modern control strategies [36–40].In the context of block-oriented models, knowledge of process dynamics can be a good starting point for identification [52]. In this regard, the Best Linear Approximation or BLA of a nonlinear system [52–55] can be used. For the specific case of Wiener-Hammerstein models, the BLA does not provide information about the dynamics of each LTI block. Therefore, all BLA-based Wiener-Hammerstein identification methods have concentrated their efforts on the BLA division to generate initial estimates and avoid suboptimal local minima in an optimization procedure. In Sjöberg et al. [45], both LTI subsystems are initialized with all possible BLA partitions and least squares optimization is applied to fit the nonlinearity. Although identification results are very good, the number of possible partitions (combinations) grows with the number of poles and zeros of the BLA and therefore the computational cost required for this method can be very high.To avoid multiple BLA divisions, in Lauwers [44] and Sjöberg et al. [45] an “advanced” method is proposed where both LTI subsystems are overparameterized with all poles and zeros of the BLA. This method is formulated as a linear-in-the-parameters total-least-squares problem for which the back LTI subsystem is inverted and basis functions are used to represent both linear subsystems. To minimize the effect of overparameterization, an order reduction technique is applied. However, since the formulation is based on neglecting the effect of disturbances, the solution is in general not consistent if there is noise on the output. In addition, the BLA is required to be invertible.Another approach to initialize Wiener-Hammerstein models is presented by Westwick and Schoukens [42], where the poles and/or zeros of the BLA are classified by using a nonlinear transformation of the input and the output residuals (quadratic/cubic BLA). On the same context of QBLA/CBLA and in line with “brute-force” method, Westwick and Schoukens [46] propose a scanning technique for a rapid evaluation of all possible BLA partitions between both LTI blocks of the Wiener-Hammerstein system. With this evaluation, the vast majority of possible partitions are discarded. Both proposals based on QBLA/CBLA show excellent results and overcome some disadvantages of “brute-force” and “advanced” methods; however, the QBLA/CBLA estimation can be difficult (high variance) since it is estimated from the BLA residuals.In a more recent work, Schoukens et al. [41] propose a more robust method based on QBLA/CBLA. Unlike the two previous proposals, the BLA is split into a nonparametric framework. This avoids mainly the parameterization of the QBLA that can be tedious given that the number of poles and zeros tends to be high. Once the front and the back dynamics of the Wiener-Hammerstein model have been identified, a parameterization of both LTI blocks is required in an additional step. This step can be complicated because a linear phase shift can be present in the nonparametric estimate.To avoid QBLA/CBLA estimation, in Vanbeylen [43] a fractional model parameterization based in multiplicities (powers) of poles and zeros is presented. The problem is formulated in the frequency domain and fractional exponents indicate which poles and zeros belong to each subsystem after an optimization problem is solved. Once the poles and zeros of the BLA have been classified, both LTI blocks must be parameterized in an additional step.All methods mentioned here, each with its advantages and disadvantages, identify Wiener-Hammerstein models from the BLA. However, all require high user interaction to parameterize the LTI blocks and/or a final optimization to refit all parameters of the Wiener-Hammerstein model. In this work, we show that it is possible to obtain a good Wiener-Hammerstein model from the BLA by solving a single optimization problem, where user interaction is only required at the beginning, just configuring simple parameters of an evolutionary algorithm.Hereinafter, this paper is organized as follows. In Section2, Wiener-Hammerstein systems formulation is revisited together with some relevant information about the BLA. The proposed evolutionary algorithm, WH-EA, is presented and described in detail in Section 3, while its application and results on a numerical example and on the benchmark data SYSID’09 are presented in Section 4. Finally, in Section 5, some conclusions are reported.


## 2. Background

### 2.1. Wiener-Hammerstein Model

A Wiener-Hammerstein model consists of two LTI subsystemsGw(z) and Gh(z) surrounding a static nonlinear function f(v(t)) (Figure 1). Both LTI subsystems can be represented in the discrete-time domain as rational transfer functions in factorized form:(1)vt=Gwzut=Kw∏i=1nbz-zwi/1-zwi∏i=1naz-pwi/1-pwiut,y^t=Ghzwt=Kh∏i=1ndz-zhi/1-zhi∏i=1ncz-phi/1-phiwt,where z is the discrete-time operator, Kw, pw1⋯pwna and zw1⋯zwnb represent the static gain, poles, and zeros of the front LTI block, respectively, and Kh, ph1⋯phnc and zh1⋯zhnd represent the static gain, poles, and zeros of the back LTI block, respectively.Figure 1
Wiener-Hammerstein model.The nonlinearity can be represented as a linear combination of a finite set(M) of basis functions:(2)wt=fvt=∑m=1Mβmfmvt,where v(t) and w(t) are the input and output of the static nonlinearity, βm are weighting parameters to be estimated, and fm are basis functions.From (1) and (2), the output of the Wiener-Hammerstein model is analytically related to the input through the following expression:(3)y^t,θ=Ghz,θhfθNL,Gwz,θwut,where(4)θw=Kw,zw1,zw2⋯zwnb,pw1,pw2⋯pwna,θNL=β1,β2,…,βm,θh=Kh,zh1,zh2⋯zhnd,ph1,ph2⋯phnc,θ=θw,θNL,θh.The challenge is to find the bestθ so that the predicted output y^(t,θ) is as close as possible to the measured output y(t). Without prior knowledge of the system, this identification problem is not easy to solve, because there are some inconveniences that must be overcome:(i)
Parametersna, nb, nc, and nd are not known (i.e., the structure of the LTI blocks is unknown).(ii)
The order and basis functions for nonlinearity are not known.(iii)
Without adequate initial values ofθ, it is quite possible that the optimization process, trying to find the best θ, gets stuck in a local minimum.(iv)
Internal variablesv(t) and w(t) are not measurable.As a complement to the formulation presented, the following assumptions are made about the system.Assumption 1.
The nonlinear system to be identified can be described by (3).Assumption 2.
The Wiener-Hammerstein system will be identified from an input/output data set{ut,yt}t=1N. The input signal u(t) is Gaussian or equivalent (see Section 2.2 for more details), while the measured output y(t) may be corrupted by stationary additive noise n(t). It is further assumed that the noise is independent of the input excitation signal:(5)yt=y0t+nt.Assumption 3.
There is no cancellation of poles and zeros and all poles of both LTI subsystems must be within the unit circle.Assumption 4.
Nonlinearity is static and its current outputwt only depends on the current input vt (i.e., the nonlinearity has no memory).


### 2.2. The Best Linear Approximation (BLA) of a Wiener-Hammerstein System

The BLA of a nonlinear system for a given class of excitation signals is a linear model that minimizes the expected mean square error between the true output of the nonlinear system and the output of the linear model [56]:(6)GBLAz=arg⁡minGz⁡Eyt-Gzut2,where ut is the input that excites the nonlinear system, y(t) is the measured output, and E is the expectation operator. An alternative way to obtain the BLA of a nonlinear system is in a nonparametric framework:(7)GBLAjωk=SyujωkSuujωk,where Syu(jωk) is the cross-power spectrum between the output y(t) and the input u(t) and Suu(jωk) is the auto power spectral density of u(t) [54, 57].The BLA depends on the excitation power spectrum (bandwidth and amplitude level) and excitation probability density function. Therefore, obtaining the BLA is restricted to the type of input signal that excites the process. Most estimation methods to obtaining the BLA use Gaussian noise signals or equivalent [58].When a nonlinear system is excited with a Gaussian signal or equivalent, according to Bussgang’s theorem [59], the nonlinearity can be replaced by a constant (KNL). Therefore, in the specific case of a Wiener-Hammerstein system, the BLA can be defined by the following expression:(8)GBLAz=KNLGwhz,where Gwh(z) represents the dynamics of the nonlinear system:(9)Gwhz=∏i=1nb+ndz-zi/1-zi∏i=1na+ncz-pi/1-piut.It is evident thatp1⋯pna+nc and z1⋯znb+nd are the poles and zeros that must be assigned to Gw(z) and Gh(z). Although the BLA does not provide information to distinguish the dynamics between both LTI subsystems, knowledge of the overall dynamics of a Wiener-Hammerstein system is a good starting point to identify such systems.


## 2.1. Wiener-Hammerstein Model

A Wiener-Hammerstein model consists of two LTI subsystemsGw(z) and Gh(z) surrounding a static nonlinear function f(v(t)) (Figure 1). Both LTI subsystems can be represented in the discrete-time domain as rational transfer functions in factorized form:(1)vt=Gwzut=Kw∏i=1nbz-zwi/1-zwi∏i=1naz-pwi/1-pwiut,y^t=Ghzwt=Kh∏i=1ndz-zhi/1-zhi∏i=1ncz-phi/1-phiwt,where z is the discrete-time operator, Kw, pw1⋯pwna and zw1⋯zwnb represent the static gain, poles, and zeros of the front LTI block, respectively, and Kh, ph1⋯phnc and zh1⋯zhnd represent the static gain, poles, and zeros of the back LTI block, respectively.Figure 1
Wiener-Hammerstein model.The nonlinearity can be represented as a linear combination of a finite set(M) of basis functions:(2)wt=fvt=∑m=1Mβmfmvt,where v(t) and w(t) are the input and output of the static nonlinearity, βm are weighting parameters to be estimated, and fm are basis functions.From (1) and (2), the output of the Wiener-Hammerstein model is analytically related to the input through the following expression:(3)y^t,θ=Ghz,θhfθNL,Gwz,θwut,where(4)θw=Kw,zw1,zw2⋯zwnb,pw1,pw2⋯pwna,θNL=β1,β2,…,βm,θh=Kh,zh1,zh2⋯zhnd,ph1,ph2⋯phnc,θ=θw,θNL,θh.The challenge is to find the bestθ so that the predicted output y^(t,θ) is as close as possible to the measured output y(t). Without prior knowledge of the system, this identification problem is not easy to solve, because there are some inconveniences that must be overcome:(i)
Parametersna, nb, nc, and nd are not known (i.e., the structure of the LTI blocks is unknown).(ii)
The order and basis functions for nonlinearity are not known.(iii)
Without adequate initial values ofθ, it is quite possible that the optimization process, trying to find the best θ, gets stuck in a local minimum.(iv)
Internal variablesv(t) and w(t) are not measurable.As a complement to the formulation presented, the following assumptions are made about the system.Assumption 1.
The nonlinear system to be identified can be described by (3).Assumption 2.
The Wiener-Hammerstein system will be identified from an input/output data set{ut,yt}t=1N. The input signal u(t) is Gaussian or equivalent (see Section 2.2 for more details), while the measured output y(t) may be corrupted by stationary additive noise n(t). It is further assumed that the noise is independent of the input excitation signal:(5)yt=y0t+nt.Assumption 3.
There is no cancellation of poles and zeros and all poles of both LTI subsystems must be within the unit circle.Assumption 4.
Nonlinearity is static and its current outputwt only depends on the current input vt (i.e., the nonlinearity has no memory).


## 2.2. The Best Linear Approximation (BLA) of a Wiener-Hammerstein System

The BLA of a nonlinear system for a given class of excitation signals is a linear model that minimizes the expected mean square error between the true output of the nonlinear system and the output of the linear model [56]:(6)GBLAz=arg⁡minGz⁡Eyt-Gzut2,where ut is the input that excites the nonlinear system, y(t) is the measured output, and E is the expectation operator. An alternative way to obtain the BLA of a nonlinear system is in a nonparametric framework:(7)GBLAjωk=SyujωkSuujωk,where Syu(jωk) is the cross-power spectrum between the output y(t) and the input u(t) and Suu(jωk) is the auto power spectral density of u(t) [54, 57].The BLA depends on the excitation power spectrum (bandwidth and amplitude level) and excitation probability density function. Therefore, obtaining the BLA is restricted to the type of input signal that excites the process. Most estimation methods to obtaining the BLA use Gaussian noise signals or equivalent [58].When a nonlinear system is excited with a Gaussian signal or equivalent, according to Bussgang’s theorem [59], the nonlinearity can be replaced by a constant (KNL). Therefore, in the specific case of a Wiener-Hammerstein system, the BLA can be defined by the following expression:(8)GBLAz=KNLGwhz,where Gwh(z) represents the dynamics of the nonlinear system:(9)Gwhz=∏i=1nb+ndz-zi/1-zi∏i=1na+ncz-pi/1-piut.It is evident thatp1⋯pna+nc and z1⋯znb+nd are the poles and zeros that must be assigned to Gw(z) and Gh(z). Although the BLA does not provide information to distinguish the dynamics between both LTI subsystems, knowledge of the overall dynamics of a Wiener-Hammerstein system is a good starting point to identify such systems.


## 3. The Evolutionary Algorithm (WH-EA)

In this paper, the identification of a Wiener-Hammerstein system is addressed as an optimization problem, which is formulated considering the following issues:(i)
The BLA is estimated in the first instance.(ii)
The poles and zeros of the BLA must be classified to find the dynamics of the front and back of the Wiener-Hammerstein model.(iii)
The pole-zero locations of the BLA can change moderately to improve modelling errors.(iv)
Without loss of generality, it is possible to model a Wiener-Hammerstein system considering that both linear blocks have unit gain (gains fromGw(z) and Gh(z) are not part of the optimization problem).(v)
The nonlinear static function is modelled as a piecewise function represented by a set of points in thev, w plane.To explain in detail how WH-EA works, this section has been divided into three parts. First part explains how Wiener-Hammerstein model is coded for individuals in the population; in addition, the optimization problem statement is presented. Second part explains in detail the customized genetic operators developed. Finally, the third part explains the general procedure of WH-EA.


### 3.1. Optimization Problem Statement

Changing pole/zero locations of the BLA to improve modelling error implies new estimates around the known values. These locations for both linear subsystems are coded in a single vector as follows:(10)P=zc1,…,zcnc,zr1,…,zrnr,zi1,…,zinc,pc1,…,pcmc,pr1,…,prmr,pi1,…,pimc,where zr1,…,zrnr and pr1,…,prmr contain the locations of the real zeros and poles, respectively, and zc1,…,zcnc and zi1,…,zinc contain the real and imaginary parts of complex conjugate zeros, respectively, while pc1,…,pcmc and pi1,…,pimc contain the real and imaginary parts of complex conjugate poles, respectively. The values of nc, nr, mc, and mr depend on the number of zeros and poles (real and/or complex conjugates) of the BLA.Poles and zeros contained in (10) must be classified to obtain the dynamics of the front and back blocks of a Wiener-Hammerstein model. This classification is performed using a binary vector:(11)C=xz1,…,xznc+nr,xp1,…,xpmc+mr.The first part of the vector,C(xz1,…,xznc+nr), is associated with zc1,…,zcnc, zr1,…,zrnr and indicates the zeros classification, while its second part, C(xp1,…,xpmc+mr), is associated with pc1,…,pcmc, pr1,…,prmr indicating the poles classification. Note that imaginary parts are not considered for classification since they are already associated with their corresponding real parts. It is assumed that if xzith=1, the corresponding ith element of P with i=1,…,nc+nr (i.e., a real zero or a pair of complex conjugated zeros) will belong to the subsystem Gw(z); otherwise, it will belong to the subsystem Gh(z). In the same way, this correspondence can be applied to classify the poles using xp1,…,xpmc+mr.For example, if a nonlinear system is approximated by a BLA with four poles,p1,2=-0.32±0.77j, p3=-0.11, and p4=0.17, and three zeros, z1,2=1.41±0.56j, z3=1.1, then nc=1, nr=1, mc=1, and mr=2, P would be structured as 1.41,1.1,0.56,-0.32,-0.11,0.17,0.77, and vector C should contain five elements whose values switch between zero and one as the algorithm evolves. By way of illustration if C=1,0,0,1,1, then Gw(z) would have two zeros and two poles: z1,2=1.41±0.56j, p3=-0.11, p4=0.17, while Gh(z) would have a zero and two poles: z3=1.1, p1,2=-0.32±0.77j.With respect to nonlinear static function, let us consider that it is represented by a set ofn points:(12)B=v1,…,vn,w1,…,wn,where the pairs (v1,w1),…,vn,wn correspond to their coordinates in a two-dimensional v-w plane. The location of these points and the interpolation method used will determine the quality of the captured static nonlinearity.The proposed evolutionary algorithm is based on stochastic population of candidate solutions (individuals). Each individual contains genetic information related to(i)
the pole/zero locations in theZ-plane of the linear subsystems (P),(ii)
the point coordinates representing the nonlinear static function (B),(iii)
and the pole/zero classification for blocksGw(z) and Gh(z) (C),such that any Wiener-Hammerstein model ((1) and (2)) can be easily described from this coded information. Recall that gains from linear blocks are assumed to be 1 and that parameters na, nb, nc, and nd will be implicitly optimized and they will depend on the structure of vector C.To find the best set of parameters, an optimization problem is stated based on a prediction-error method and the typical mean-squared error criterion (although any other criteria can be used in the proposed method, such as the mean absolute or maximum error criteria):(13)εt,θ=yt-y^t,θ,(14)Jθ=1N∑t=1Nε2t,θ,where θ=P,B,C and the solution of the optimization problem is stated as(15)θ^=argminθJθ,where θ^ contains the genetic information from the best individual at the end of generations.


### 3.2. Genetic Operators

Customized mutation and crossover operators will be developed taking in mind the problem at hand: to identify all parameters of the Wiener-Hammerstein model in a single optimization trial. Figure2 shows the structure of an individual as well as the genetic operators developed on each piece of genetic information. Note that i and g have been introduced into the formulation. Subscript i represents an individual in the population, while the superscript g indicates the current population.Figure 2
Structure of individual and genetic operations performed on each piece of genetic information.The specific mutation and crossover operators designed are randomly selected to maintain a balance between exploration and exploitation of the search space. Mutation operations are used to maintain genetic diversity, while crossover operations allow genetic information from the best individuals to be combined and disseminated throughout the generations. Further details on how the algorithm works will be given in Section3.3.


#### 3.2.1. Location in the Z-Plane of Poles and Zeros

Theoretically in a Wiener-Hammerstein model, the pole-zero locations ofGw(z) and Gh(z) subsystems correspond to the pole-zero locations of the BLA; however, it is well known that once the BLA has been divided, a refit can be used to improve the modelling error. In this regard, the proposed algorithm considers that while the BLA is divided and nonlinearity is captured, the pole-zero locations can change subtly.Both operations used on this portion of genetic information produce offspring vectorP~g, which directly inherits from its parent Pig all the genetic information except in a gene. This gene will be selected using a random integer number rzp∈1,…,nr+2nc+mr+2mc and modified according to the corresponding genetic operatormutation M.1 orcrossover C.1.Mutation M.1. The selected gene is mutated to explore in an individualized way new pole-zero locations of the BLA. A new location P~jg is determined by a random number Nzp with Gaussian distribution:(16)P~jg=Pi,jg+Nzp0,σ2gifj=rzpPi,jgotherwise,where j=1⋯nr+2nc+mr+2mc. Pi,jg and P~jg represent thejth elements of vectors Pig and P~g, respectively.The new locations for poles and zeros are explored within a search space defined byPmin and Pmax; therefore, Pjmin≤P~jg≤Pjmax, where Pjmin and Pjmax are thejth elements of vectors Pmin and Pmax, respectively (see Section 3.3 for more details on search space for poles and zeros).Aggressiveness of mutations can be controlled through the standard deviation:(17)σ2g=Δs100σini21+g∗σratio2,σratio2=σini2/σend22-1MaxGen-1,where MaxGen is the predefined number of algorithm generations; σratio2 is the rate at which the standard deviation will decrease from σini2 to σend2 as the generations pass (see Figure 3); the parameter Δs bounds the limits of the interval in which the selected gene can be moved. For this mutation, Δs=Przpmax-Przpmin. Variation of σ2(g) will allow mutations to be more subtle in the last generations to achieve a fine-tuning of the corresponding parameters.Figure 3
Variation of standard deviation over generations to control the aggressiveness of mutations.Crossover C.1. The selected gene is formed using genetic information from the parent, Pi,jg, combined with the corresponding genetic information from the best individual, Pbest,jg, in the current population:(18)P~jg=Pi,jg+Pbest,jg2ifj=rzpPi,jgotherwise,where j=1⋯nr+2nc+mr+2mc.


#### 3.2.2. Nonlinear Static Function

As the algorithm evolves, points for nonlinear static function must be located adequately in thev-w plane. Here any type of interpolation can be used to capture the static nonlinearity. To achieve a good fit, mutations M.2 and M.3 plus a crossover operation are used. Both mutations used on this portion of genetic information produce offspring vector B~g, which directly inherits from its parent Big all the genetic information except in two genes. This pair of genes represents the coordinates of the point that will be modified. Unlike mutation operations, the crossover operation generates offspring with a single modified gene which corresponds to the ordinate of a point. Given the correspondence between the abscissa and the ordinate of a point, for the three operations a single integer random number rbp∈1,n will allow us to select the gene(s) to be modified.Mutation M.2. This genetic operation allows us to explore in the v-w plane new positions for the points. The mutation in both genes is handled by random numbers (Nv,Nw) with Gaussian distribution:(19)B~jg=Bi,jg+Nv0,σ2gifj=rbpBi,jg+Nw0,σ2gifj=n+rbpBi,jgotherwisewith j=1⋯2n. Bi,jg and B~jg represent the jth elements of vectors Big and B~g, respectively. To avoid overlapping points, bounds for mutations on the abscissa axis are set depending on the selected point to mutate (B~rbpg) and the location of its neighbors according to(i)
Bi,rbpg+α<B~rbpg<Bi,rbp+1g-α; if rbp=1,(ii)
Bi,rbp-1g+α<B~rbpg<Bi,rbp+1g-α; if rbp=2⋯n-1,(iii)
Bi,rbp-1g+α<B~rbpg<Bi,rbpg-α; if rbp=n,whereα is a user-defined parameter that indicates how close the points can be located. To achieve a good fit of the nonlinearity α must be small, relative to the search space on the abscissa axis defined by vmin and vmax. Note that the horizontal boundaries for the endpoints are delimited by their position and their left or right neighbor, respectively. Bounds for mutations on ordinate axis are fixed and equal for all points. This allows each point to move freely throughout the search space on the ordinate axis defined by wmin and wmax. The vertical and horizontal bounds for a selected point are illustrated in Figure 4. When mutation M.2 is required, the selected point can be changed to a new position within the grey rectangle. Details on how to determine vmin, vmaxwmin, and wmax will be given in Section 3.3.Figure 4
Bounds for mutation M.2. Grey area indicatesB~rbpg and B~n+rbpg feasible space.To achieve a fine-tuning of all nonlinearity points, mutations’ aggressiveness can be controlled through standard deviation (17) as in mutation M.1. Note that Δs=wmax-wmin is constant for all mutations over ordinate axis, while for abscissa axis mutations, Δs can be calculated as(20)Δs=Bi,rbp+1g-α-vmin;ifrbp=1Bi,rbp+1g-Bi,rbp-1g-2αifrbp=2⋯n-1vmax-Bi,rbp-1g-αifrbp=n.Mutation M.2 has great potential to explore the searching space. This genetic operation will locate points where there are slope changes. During first generations, it is useful to shape nonlinearity, while in last ones, it allows a refinement. However, when a point is located where there is a slope change, it could be kept in this location until the end of generations, especially when there are abrupt changes in the slope. Because jumps between points are not allowed with mutation M.2, when a point is kept in a place where there is a significant change of slope, one or more points would remain trapped to the left or right of it. This would lead to having redundant points in a segment that would not require so many or worse, to having a segment (curvature) that would not contain enough points. To avoid this drawback, the exploration in the search space is complemented with mutation M.3.Mutation M.3. This genetic operation is designed to concentrate as many points as possible on the curvatures that nonlinearity can have. Therefore, it will be required that each point can be displaced on the abscissa axis by jumping one or more positions of the other points. Let us define a segment as the horizontal space between two consecutive points (so for n points there will be n-1 segments); then a random integer number rs∈1,n-1 will indicate to which segment the selected point will move. The first half of the offspring vector is found using the following expression:(21)B~jg=Bi,rsg+Bi,rs+1g2ifj=rbpBi,jgotherwisewith j=1⋯n. Note that if rbp=rs or rbp=rs-1, the corresponding point will not make a jump but it will be located at the midpoint between its current position and the position of the point on the right or left, respectively. As in the mutation M.2 to prevent points from getting too close together, the α parameter is also used in this mutation; therefore, a jump is conditioned to the space available in the selected segment to accommodate a new point. Minimum space should be 2α. If this condition is not met, rs must be regenerated to randomly search for another segment.To provide a smooth transition between adjacent segments, gene mutation corresponding to the position on the ordinates axis is performed using a quadratic interpolation. To do that, three neighboring points are required. The second half of the offspring vector is found using the following expression:(22)B~jg=B~j-ng2,B~j-ng,1∗k2k1k0ifj=n+rbpBi,jgotherwisewith j=n+1⋯2n. B~j-ng is the v-coordinate of the selected point to mute which can be found with (21); k0, k1, and k2 are the coefficients of the quadratic polynomial Ψ defined by three adjacent points selected once the new v-coordinate of the point that is mutating is known. The three adjacent points can be selected directly when a point has mutated to the first or last segment,(23)Ψ=fBi,1g,Bi,n+1g;Bi,2g,Bi,n+2g;Bi,3g,Bi,n+3g:ifrs=1fBi,ng,Bi,2ng;Bi,n-1g,Bi,2n-1g;Bi,n-2g,Bi,2n-2g:ifrs=n-1,while if the point has mutated to a nonextreme segment, the three adjacent points can be selected using the two points that define that segment plus one on its right or left. For more effective exploration, a random number r3p∈(0,1] is used for selection:(24)Ψ=fBi,rs-1g,Bi,n+rs-1g;Bi,rsg,Bi,n+rsg;Bi,rs+1g,Bi,n+rs+1g:ifr3p≤0.5fBi,rsg,Bi,n+rsg;Bi,rs+1g,Bi,n+rs+1g;Bi,rs+2g,Bi,n+rs+2g:otherwise.Figure5 illustrates how a jump occurs with mutation M.3. Notice that the new ordinate is calculated according to the polynomial formed by the two points of the segment plus a point to the right: that is, r3p>0.5. After a jump has occurred, an ascending reordering of the points with respect to the abscissa values is necessary.Figure 5
Mutation M.3 withrbp=2, rs=4, and r3p>0.5. Jump to the selected segment (dashed line). Quadratic polynomial (solid line).Crossover C.2. This genetic operation works just like crossover C.1 and is applied only to vary the position of a point on the ordinate axis:(25)B~jg=Bi,jg+Bbest,jg2ifj=n+rbpBi,jgotherwisewith j=n+1⋯2n. Bbest,jg is the jth element of vector Bbestg, which corresponds to the individual in the current population g with the best fitness value.


#### 3.2.3. Pole-Zero Classification

Due to the stochastic nature of evolutionary algorithms, the binary values of (11) will change as the algorithm evolves, generating different structures of Gw(z) and Gh(z). The evolution of this piece of genetic information is handled by a simple mutation operator.Mutation M.4. Unlike the previous ones, this operator generates a new vector C~g that depends entirely on the effects of mutation, meaning that for this piece of genetic information there is no information exchange between generations. This allows free testing of different structures for Gw(z) and Gh(z) to avoid premature convergence. When this operation is required, a random process will generate the mutation vector:(26)C~jg=1ifNc≤0.50otherwisewith j=1⋯nc+nr+mc+mr. C~jg is the jth element of vector C~g. Nc is a random number with standard uniform distribution on the open interval (0,1). Note that the structure of C~g is built under two considerations: the LTI subsystems cannot be improper and the sum of zeros and the sum of the poles between both subsystems must be equal to the number of zeros and poles of the BLA, respectively.


### 3.3. WH-EA Description

WH-EA is an elitist evolutionary algorithm that evolves a population ofNP individuals. Each individual contains three portions of genetic information related to the parameters of a Wiener-Hammerstein model θig=Pig,Big,Cig. Like any other evolutionary algorithm, WH-EA is inspired by biological evolution over generations. Starting from an initial population, new generations are created using information of the current generation g and performing crossover and/or mutation operations and selection based on the fitness of the new individuals. Algorithm 1 shows a pseudocode of main steps performed in WH-EA, whereas details of all the parameters that the algorithm uses are shown in WH-EA Parameters.Algorithm 1:Pseudocode of WH-EA.
(1) Initialise the population;
(2) Evaluate fitness of all population;
(3) for  g=1  toMaxGendo
(4)Find  θbestg
(5)Random selection of a individual (r1);
(6)Compute γ(g);
(7)if  rpznl≤γ(g)   then
(8)Compute B~g using Algorithm 2;
(9)else
(10)Compute P~g using Algorithm 3;
(11)end if
(12)if  rc≤ξthen
(13)Compute C~g using Mutation M.4;
(14)end if
(15)Update population;
(16) end for
(17) Print θbestMaxGenAlgorithm 2:Modify two-dimensional points for nonlinear function.
(1)if  rnmc≤δnlthen
(2)Compute η(g)
(3)if  rmm≤ηmin+η(g)then
(4)Mutation M.2;
(5)else
(6)Mutation M.3;
(7)end if
(8)else
(9)Crossover C.2;
(10)end ifAlgorithm 3:Modify polo/zero locations.
(1)if  rlmc≤δzpthen
(2)Mutation M.1;
(3)else
(4)Crossover C.1;
(5)end ifInitialize the Population. The initial population P0B0C0 contains NP individuals generated within a search space delimited by lower and upper bounds. For each piece of genetic information, the lower and upper bounds can be determined from the information provided by the BLA (location of poles and zeros and static gain).Real and imaginary values of poles and zeros from the BLA without modifications are introduced as part of the first individualP10 of the initial population. For the rest, mutation M.1 (16) is used but fixing the parent vector as P10, that is, the BLA.Therefore, mutation M.1 must be executedNP-1 times under the above conditions to generate mutated versions of the BLA. Initialization and mutation M.1 are performed considering a search space delimited by lower bound Pmin and upper bound Pmax which are set around the values of the BLA:(27)Pjmin=P1,j0-Υjmin,Pjmax=P1,j0+Υjmax,where Υjmin and Υjmax represent the jth elements of the user-defined vectors Υmin and Υmax, respectively, indicating how much the BLA’s pole/zero position can change as the algorithm evolves.On the other hand, the initial population corresponding to two-dimensional points must be generated within a search space defined horizontally by the minimum (vmin) and maximum (vmax) amplitude of v(t) (the output of Gw(z)) and vertically by the minimum (wmin) and maximum (wmax) amplitude of w(t) (the input to Gh(z)). Although v(t) and w(t) are not known, vmin and vmax depend on the input u(t) and Gw(z). Since Gw(z) is a linear system, the following expressions can be used:(28)vmin=Ω∗umin,vmax=Ω∗umax,where umin and umax are the minimum and maximum values of the signal u(t), respectively, and Ω is a scaling factor which depends on the pole/zero locations and the static gain of Gw(z). Without loss of generality, it is possible to model a Wiener-Hammerstein system considering that both linear blocks have unit gain.Regardingwmin and wmax, if the input u(t) with mean umean enters Gw(z) regardless of its structure, the mean of the output signal vmean will be equal to umean. This is not the case when the signal v(t) enters the nonlinear block since the gain KNL might provide an offset to w(t). However, since Gh(z) is a linear block, the mean of the output ymean will be equal to the mean of the incoming signal wmean. With this information, the straight line of Figure 6 can be drawn, and wmin and wmax can be found using the following equations:(29)wmin=ymean+KNLΩumin-umean,wmax=ymean+KNLΩumax-umean.Figure 6
Information to define the search space for the nonlinear function.As can be seen from (28) to (29), the search space for nonlinearity depends on the input signal, output signal, static gain of the BLA, and Ω which is a user-defined parameter. Since both linear subsystems will be estimated with unit gain, neither of these will amplify their input signals, therefore vmin>umin and vmax<umax. For these two conditions to be met, Ω must be less than one. In the same way, it must be observed that wmin<ymin and wmax>ymax. If Ω is less than one, the first pair of conditions will always be met; however, there is no guarantee that the second pair of conditions will be met. Since it is possible to perform this check prior to the execution of the algorithm, if the second pair of conditions are not met, Ω must be increased, but considering that it must be less than one. It should also be taken into account that if Ω is too large, the search space will be larger than necessary, so the algorithm will cost more to estimate static nonlinearity.To initialize this portion of genetic information, then points are uniformly distributed between Ωumin and Ωumax and located on the straight line shown in Figure 6. These points are introduced as part of the first individual in the population B10. The corresponding genetic information for the rest of individuals is generated using mutation M.2 (19), but considering that the parent vector is always B10.Finally, genetic information corresponding to pole-zero classification is initialized directly using mutation M.4 (26) NP times.Evaluate Fitness. Performance of each individual in the population is defined by a fitness criterion which can be calculated using (14), where θ is obtained from the encoded information in Pig,Big,Cig.The Offspring. Once population has been initialized, for each generation a random integer number r1∈1,NP will be used to select the parent from which an offspring P~g will be generated. As can be seen in Algorithm 1, not all genetic operators are applied at the same time to generate offspring; this can help to expand diversity and avoid premature convergence.One or two pieces of the offspring genetic information will be randomly selected for modification according to their respective genetic operators. A random numberrpznl∈(0,1] chooses between modifying the portion related to static nonlinearity using Algorithm 2 or the portion of genetic information related to pole/zero locations using Algorithm 3. The probability for this selection is handled by the control parameter γ(g) defined as(30)γg=γini1+g∗γrat,γrat=γini/γend2-1MaxGen-1,where γrat is the rate at which the probability γ(g) will decrease from initial probability γini to final probability γend as generations pass; therefore, 0<γend<γini≤1. If γini=1, the probability of modifying the genetic information of nonlinearity in the first generations will be high, while the probability of modifying the location of poles and zeros will be low. On the other hand, if and γend=0.5, in the final generations, the algorithm will modify with equal probability both portions of genetic information. The selection of these values is justified by the fact that pole/zero locations are known and they will only be fine-tuned within a suitable search space to amend possible errors in the BLA estimation, whereas nonlinearity is completely unknown, so the algorithm should focus more on this portion of genetic information during first generations.Variation of genetic information corresponding to the classification of poles and zeros for both LTI subsystems is handled by a comparison between a random numberrc∈(0,1] and the probability ξ∈(0,1]. The value of probability ξ is defined by the user and will be constant throughout the evolution of the algorithm. Figure 7 shows the behaviour of the control parameters (probabilities) used to select the portions of genetic information that will be modified in each generation.Figure 7
Control parameters of WH-EA for selection of the genetic information to be modified in generationg.Algorithm2 is used to modify the genetic information related to nonlinear static function. The control parameter δnl∈(0,1] indicates the probability with which the mutation (either M.2 or M.3) or crossover C.2 will be used. Probability of selecting M.2 or M.3 is variable with respect to the generations. During first generations, mutation M.3 is not necessary, since the nonlinearity can be captured thanks to the two-dimensional points movements due to mutation M.2 and crossover C.2 operations. Since it is very likely that nonlinearity includes one or more curvatures, as the algorithm evolves mutation M.3 will be required to concentrate as many points as possible on these curvatures. The variable probability for selection between both mutations is defined by(31)ηg=1-ηmin-1-ηmingMaxGen,where ηmin∈(0,0.5] is a user-defined parameter indicating the minimum probability with which the mutation M.2 can be selected. Note that according to (31) and Algorithm 2, the maximum probability is 1 and occurs in the first generation. As the algorithm evolves, this probability will decrease linearly until it reaches ηmin in the last generation. When Algorithm 2 is required a random number rnmc∈(0,1] will allow us to select either a mutation or crossover C.2. If a mutation is selected, a new random number rmm∈(0,1] will allow us to select between mutation M.2 or mutation M.3.On the other hand, Algorithm3 is used to modify the genetic information related to pole/zero locations using mutation M.1 or crossover C.1. The control parameter δzp∈(0,1] indicates the probability with which each genetic operation will be used. Since crossover C.1 causes offspring to inherit genetic information from the best individual, a small value of δzp may lead to premature convergence, whereas a value closer to 1 will cause the algorithm to converge very slowly. When Algorithm 3 is required, a random number rlmc∈(0,1] will determine the genetic operation to be used.Update. It is based on a competition between the generated offspring and the individuals of the population. The contestant with the best fitness will be the one who wins the competition. From a randomly selected individual, the offspring starts to compete until defeating an individual; when this happens the descendant will take his place in the population and the algorithm continues with the next generation. If the offspring comes to compete with all individuals and could not win, this will be discarded and the algorithm will pass to the next generation.


## 3.1. Optimization Problem Statement

Changing pole/zero locations of the BLA to improve modelling error implies new estimates around the known values. These locations for both linear subsystems are coded in a single vector as follows:(10)P=zc1,…,zcnc,zr1,…,zrnr,zi1,…,zinc,pc1,…,pcmc,pr1,…,prmr,pi1,…,pimc,where zr1,…,zrnr and pr1,…,prmr contain the locations of the real zeros and poles, respectively, and zc1,…,zcnc and zi1,…,zinc contain the real and imaginary parts of complex conjugate zeros, respectively, while pc1,…,pcmc and pi1,…,pimc contain the real and imaginary parts of complex conjugate poles, respectively. The values of nc, nr, mc, and mr depend on the number of zeros and poles (real and/or complex conjugates) of the BLA.Poles and zeros contained in (10) must be classified to obtain the dynamics of the front and back blocks of a Wiener-Hammerstein model. This classification is performed using a binary vector:(11)C=xz1,…,xznc+nr,xp1,…,xpmc+mr.The first part of the vector,C(xz1,…,xznc+nr), is associated with zc1,…,zcnc, zr1,…,zrnr and indicates the zeros classification, while its second part, C(xp1,…,xpmc+mr), is associated with pc1,…,pcmc, pr1,…,prmr indicating the poles classification. Note that imaginary parts are not considered for classification since they are already associated with their corresponding real parts. It is assumed that if xzith=1, the corresponding ith element of P with i=1,…,nc+nr (i.e., a real zero or a pair of complex conjugated zeros) will belong to the subsystem Gw(z); otherwise, it will belong to the subsystem Gh(z). In the same way, this correspondence can be applied to classify the poles using xp1,…,xpmc+mr.For example, if a nonlinear system is approximated by a BLA with four poles,p1,2=-0.32±0.77j, p3=-0.11, and p4=0.17, and three zeros, z1,2=1.41±0.56j, z3=1.1, then nc=1, nr=1, mc=1, and mr=2, P would be structured as 1.41,1.1,0.56,-0.32,-0.11,0.17,0.77, and vector C should contain five elements whose values switch between zero and one as the algorithm evolves. By way of illustration if C=1,0,0,1,1, then Gw(z) would have two zeros and two poles: z1,2=1.41±0.56j, p3=-0.11, p4=0.17, while Gh(z) would have a zero and two poles: z3=1.1, p1,2=-0.32±0.77j.With respect to nonlinear static function, let us consider that it is represented by a set ofn points:(12)B=v1,…,vn,w1,…,wn,where the pairs (v1,w1),…,vn,wn correspond to their coordinates in a two-dimensional v-w plane. The location of these points and the interpolation method used will determine the quality of the captured static nonlinearity.The proposed evolutionary algorithm is based on stochastic population of candidate solutions (individuals). Each individual contains genetic information related to(i)
the pole/zero locations in theZ-plane of the linear subsystems (P),(ii)
the point coordinates representing the nonlinear static function (B),(iii)
and the pole/zero classification for blocksGw(z) and Gh(z) (C),such that any Wiener-Hammerstein model ((1) and (2)) can be easily described from this coded information. Recall that gains from linear blocks are assumed to be 1 and that parameters na, nb, nc, and nd will be implicitly optimized and they will depend on the structure of vector C.To find the best set of parameters, an optimization problem is stated based on a prediction-error method and the typical mean-squared error criterion (although any other criteria can be used in the proposed method, such as the mean absolute or maximum error criteria):(13)εt,θ=yt-y^t,θ,(14)Jθ=1N∑t=1Nε2t,θ,where θ=P,B,C and the solution of the optimization problem is stated as(15)θ^=argminθJθ,where θ^ contains the genetic information from the best individual at the end of generations.


## 3.2. Genetic Operators

Customized mutation and crossover operators will be developed taking in mind the problem at hand: to identify all parameters of the Wiener-Hammerstein model in a single optimization trial. Figure2 shows the structure of an individual as well as the genetic operators developed on each piece of genetic information. Note that i and g have been introduced into the formulation. Subscript i represents an individual in the population, while the superscript g indicates the current population.Figure 2
Structure of individual and genetic operations performed on each piece of genetic information.The specific mutation and crossover operators designed are randomly selected to maintain a balance between exploration and exploitation of the search space. Mutation operations are used to maintain genetic diversity, while crossover operations allow genetic information from the best individuals to be combined and disseminated throughout the generations. Further details on how the algorithm works will be given in Section3.3.


### 3.2.1. Location in the Z-Plane of Poles and Zeros

Theoretically in a Wiener-Hammerstein model, the pole-zero locations ofGw(z) and Gh(z) subsystems correspond to the pole-zero locations of the BLA; however, it is well known that once the BLA has been divided, a refit can be used to improve the modelling error. In this regard, the proposed algorithm considers that while the BLA is divided and nonlinearity is captured, the pole-zero locations can change subtly.Both operations used on this portion of genetic information produce offspring vectorP~g, which directly inherits from its parent Pig all the genetic information except in a gene. This gene will be selected using a random integer number rzp∈1,…,nr+2nc+mr+2mc and modified according to the corresponding genetic operatormutation M.1 orcrossover C.1.Mutation M.1. The selected gene is mutated to explore in an individualized way new pole-zero locations of the BLA. A new location P~jg is determined by a random number Nzp with Gaussian distribution:(16)P~jg=Pi,jg+Nzp0,σ2gifj=rzpPi,jgotherwise,where j=1⋯nr+2nc+mr+2mc. Pi,jg and P~jg represent thejth elements of vectors Pig and P~g, respectively.The new locations for poles and zeros are explored within a search space defined byPmin and Pmax; therefore, Pjmin≤P~jg≤Pjmax, where Pjmin and Pjmax are thejth elements of vectors Pmin and Pmax, respectively (see Section 3.3 for more details on search space for poles and zeros).Aggressiveness of mutations can be controlled through the standard deviation:(17)σ2g=Δs100σini21+g∗σratio2,σratio2=σini2/σend22-1MaxGen-1,where MaxGen is the predefined number of algorithm generations; σratio2 is the rate at which the standard deviation will decrease from σini2 to σend2 as the generations pass (see Figure 3); the parameter Δs bounds the limits of the interval in which the selected gene can be moved. For this mutation, Δs=Przpmax-Przpmin. Variation of σ2(g) will allow mutations to be more subtle in the last generations to achieve a fine-tuning of the corresponding parameters.Figure 3
Variation of standard deviation over generations to control the aggressiveness of mutations.Crossover C.1. The selected gene is formed using genetic information from the parent, Pi,jg, combined with the corresponding genetic information from the best individual, Pbest,jg, in the current population:(18)P~jg=Pi,jg+Pbest,jg2ifj=rzpPi,jgotherwise,where j=1⋯nr+2nc+mr+2mc.


### 3.2.2. Nonlinear Static Function

As the algorithm evolves, points for nonlinear static function must be located adequately in thev-w plane. Here any type of interpolation can be used to capture the static nonlinearity. To achieve a good fit, mutations M.2 and M.3 plus a crossover operation are used. Both mutations used on this portion of genetic information produce offspring vector B~g, which directly inherits from its parent Big all the genetic information except in two genes. This pair of genes represents the coordinates of the point that will be modified. Unlike mutation operations, the crossover operation generates offspring with a single modified gene which corresponds to the ordinate of a point. Given the correspondence between the abscissa and the ordinate of a point, for the three operations a single integer random number rbp∈1,n will allow us to select the gene(s) to be modified.Mutation M.2. This genetic operation allows us to explore in the v-w plane new positions for the points. The mutation in both genes is handled by random numbers (Nv,Nw) with Gaussian distribution:(19)B~jg=Bi,jg+Nv0,σ2gifj=rbpBi,jg+Nw0,σ2gifj=n+rbpBi,jgotherwisewith j=1⋯2n. Bi,jg and B~jg represent the jth elements of vectors Big and B~g, respectively. To avoid overlapping points, bounds for mutations on the abscissa axis are set depending on the selected point to mutate (B~rbpg) and the location of its neighbors according to(i)
Bi,rbpg+α<B~rbpg<Bi,rbp+1g-α; if rbp=1,(ii)
Bi,rbp-1g+α<B~rbpg<Bi,rbp+1g-α; if rbp=2⋯n-1,(iii)
Bi,rbp-1g+α<B~rbpg<Bi,rbpg-α; if rbp=n,whereα is a user-defined parameter that indicates how close the points can be located. To achieve a good fit of the nonlinearity α must be small, relative to the search space on the abscissa axis defined by vmin and vmax. Note that the horizontal boundaries for the endpoints are delimited by their position and their left or right neighbor, respectively. Bounds for mutations on ordinate axis are fixed and equal for all points. This allows each point to move freely throughout the search space on the ordinate axis defined by wmin and wmax. The vertical and horizontal bounds for a selected point are illustrated in Figure 4. When mutation M.2 is required, the selected point can be changed to a new position within the grey rectangle. Details on how to determine vmin, vmaxwmin, and wmax will be given in Section 3.3.Figure 4
Bounds for mutation M.2. Grey area indicatesB~rbpg and B~n+rbpg feasible space.To achieve a fine-tuning of all nonlinearity points, mutations’ aggressiveness can be controlled through standard deviation (17) as in mutation M.1. Note that Δs=wmax-wmin is constant for all mutations over ordinate axis, while for abscissa axis mutations, Δs can be calculated as(20)Δs=Bi,rbp+1g-α-vmin;ifrbp=1Bi,rbp+1g-Bi,rbp-1g-2αifrbp=2⋯n-1vmax-Bi,rbp-1g-αifrbp=n.Mutation M.2 has great potential to explore the searching space. This genetic operation will locate points where there are slope changes. During first generations, it is useful to shape nonlinearity, while in last ones, it allows a refinement. However, when a point is located where there is a slope change, it could be kept in this location until the end of generations, especially when there are abrupt changes in the slope. Because jumps between points are not allowed with mutation M.2, when a point is kept in a place where there is a significant change of slope, one or more points would remain trapped to the left or right of it. This would lead to having redundant points in a segment that would not require so many or worse, to having a segment (curvature) that would not contain enough points. To avoid this drawback, the exploration in the search space is complemented with mutation M.3.Mutation M.3. This genetic operation is designed to concentrate as many points as possible on the curvatures that nonlinearity can have. Therefore, it will be required that each point can be displaced on the abscissa axis by jumping one or more positions of the other points. Let us define a segment as the horizontal space between two consecutive points (so for n points there will be n-1 segments); then a random integer number rs∈1,n-1 will indicate to which segment the selected point will move. The first half of the offspring vector is found using the following expression:(21)B~jg=Bi,rsg+Bi,rs+1g2ifj=rbpBi,jgotherwisewith j=1⋯n. Note that if rbp=rs or rbp=rs-1, the corresponding point will not make a jump but it will be located at the midpoint between its current position and the position of the point on the right or left, respectively. As in the mutation M.2 to prevent points from getting too close together, the α parameter is also used in this mutation; therefore, a jump is conditioned to the space available in the selected segment to accommodate a new point. Minimum space should be 2α. If this condition is not met, rs must be regenerated to randomly search for another segment.To provide a smooth transition between adjacent segments, gene mutation corresponding to the position on the ordinates axis is performed using a quadratic interpolation. To do that, three neighboring points are required. The second half of the offspring vector is found using the following expression:(22)B~jg=B~j-ng2,B~j-ng,1∗k2k1k0ifj=n+rbpBi,jgotherwisewith j=n+1⋯2n. B~j-ng is the v-coordinate of the selected point to mute which can be found with (21); k0, k1, and k2 are the coefficients of the quadratic polynomial Ψ defined by three adjacent points selected once the new v-coordinate of the point that is mutating is known. The three adjacent points can be selected directly when a point has mutated to the first or last segment,(23)Ψ=fBi,1g,Bi,n+1g;Bi,2g,Bi,n+2g;Bi,3g,Bi,n+3g:ifrs=1fBi,ng,Bi,2ng;Bi,n-1g,Bi,2n-1g;Bi,n-2g,Bi,2n-2g:ifrs=n-1,while if the point has mutated to a nonextreme segment, the three adjacent points can be selected using the two points that define that segment plus one on its right or left. For more effective exploration, a random number r3p∈(0,1] is used for selection:(24)Ψ=fBi,rs-1g,Bi,n+rs-1g;Bi,rsg,Bi,n+rsg;Bi,rs+1g,Bi,n+rs+1g:ifr3p≤0.5fBi,rsg,Bi,n+rsg;Bi,rs+1g,Bi,n+rs+1g;Bi,rs+2g,Bi,n+rs+2g:otherwise.Figure5 illustrates how a jump occurs with mutation M.3. Notice that the new ordinate is calculated according to the polynomial formed by the two points of the segment plus a point to the right: that is, r3p>0.5. After a jump has occurred, an ascending reordering of the points with respect to the abscissa values is necessary.Figure 5
Mutation M.3 withrbp=2, rs=4, and r3p>0.5. Jump to the selected segment (dashed line). Quadratic polynomial (solid line).Crossover C.2. This genetic operation works just like crossover C.1 and is applied only to vary the position of a point on the ordinate axis:(25)B~jg=Bi,jg+Bbest,jg2ifj=n+rbpBi,jgotherwisewith j=n+1⋯2n. Bbest,jg is the jth element of vector Bbestg, which corresponds to the individual in the current population g with the best fitness value.


### 3.2.3. Pole-Zero Classification

Due to the stochastic nature of evolutionary algorithms, the binary values of (11) will change as the algorithm evolves, generating different structures of Gw(z) and Gh(z). The evolution of this piece of genetic information is handled by a simple mutation operator.Mutation M.4. Unlike the previous ones, this operator generates a new vector C~g that depends entirely on the effects of mutation, meaning that for this piece of genetic information there is no information exchange between generations. This allows free testing of different structures for Gw(z) and Gh(z) to avoid premature convergence. When this operation is required, a random process will generate the mutation vector:(26)C~jg=1ifNc≤0.50otherwisewith j=1⋯nc+nr+mc+mr. C~jg is the jth element of vector C~g. Nc is a random number with standard uniform distribution on the open interval (0,1). Note that the structure of C~g is built under two considerations: the LTI subsystems cannot be improper and the sum of zeros and the sum of the poles between both subsystems must be equal to the number of zeros and poles of the BLA, respectively.


## 3.2.1. Location in the Z-Plane of Poles and Zeros

Theoretically in a Wiener-Hammerstein model, the pole-zero locations ofGw(z) and Gh(z) subsystems correspond to the pole-zero locations of the BLA; however, it is well known that once the BLA has been divided, a refit can be used to improve the modelling error. In this regard, the proposed algorithm considers that while the BLA is divided and nonlinearity is captured, the pole-zero locations can change subtly.Both operations used on this portion of genetic information produce offspring vectorP~g, which directly inherits from its parent Pig all the genetic information except in a gene. This gene will be selected using a random integer number rzp∈1,…,nr+2nc+mr+2mc and modified according to the corresponding genetic operatormutation M.1 orcrossover C.1.Mutation M.1. The selected gene is mutated to explore in an individualized way new pole-zero locations of the BLA. A new location P~jg is determined by a random number Nzp with Gaussian distribution:(16)P~jg=Pi,jg+Nzp0,σ2gifj=rzpPi,jgotherwise,where j=1⋯nr+2nc+mr+2mc. Pi,jg and P~jg represent thejth elements of vectors Pig and P~g, respectively.The new locations for poles and zeros are explored within a search space defined byPmin and Pmax; therefore, Pjmin≤P~jg≤Pjmax, where Pjmin and Pjmax are thejth elements of vectors Pmin and Pmax, respectively (see Section 3.3 for more details on search space for poles and zeros).Aggressiveness of mutations can be controlled through the standard deviation:(17)σ2g=Δs100σini21+g∗σratio2,σratio2=σini2/σend22-1MaxGen-1,where MaxGen is the predefined number of algorithm generations; σratio2 is the rate at which the standard deviation will decrease from σini2 to σend2 as the generations pass (see Figure 3); the parameter Δs bounds the limits of the interval in which the selected gene can be moved. For this mutation, Δs=Przpmax-Przpmin. Variation of σ2(g) will allow mutations to be more subtle in the last generations to achieve a fine-tuning of the corresponding parameters.Figure 3
Variation of standard deviation over generations to control the aggressiveness of mutations.Crossover C.1. The selected gene is formed using genetic information from the parent, Pi,jg, combined with the corresponding genetic information from the best individual, Pbest,jg, in the current population:(18)P~jg=Pi,jg+Pbest,jg2ifj=rzpPi,jgotherwise,where j=1⋯nr+2nc+mr+2mc.


## 3.2.2. Nonlinear Static Function

As the algorithm evolves, points for nonlinear static function must be located adequately in thev-w plane. Here any type of interpolation can be used to capture the static nonlinearity. To achieve a good fit, mutations M.2 and M.3 plus a crossover operation are used. Both mutations used on this portion of genetic information produce offspring vector B~g, which directly inherits from its parent Big all the genetic information except in two genes. This pair of genes represents the coordinates of the point that will be modified. Unlike mutation operations, the crossover operation generates offspring with a single modified gene which corresponds to the ordinate of a point. Given the correspondence between the abscissa and the ordinate of a point, for the three operations a single integer random number rbp∈1,n will allow us to select the gene(s) to be modified.Mutation M.2. This genetic operation allows us to explore in the v-w plane new positions for the points. The mutation in both genes is handled by random numbers (Nv,Nw) with Gaussian distribution:(19)B~jg=Bi,jg+Nv0,σ2gifj=rbpBi,jg+Nw0,σ2gifj=n+rbpBi,jgotherwisewith j=1⋯2n. Bi,jg and B~jg represent the jth elements of vectors Big and B~g, respectively. To avoid overlapping points, bounds for mutations on the abscissa axis are set depending on the selected point to mutate (B~rbpg) and the location of its neighbors according to(i)
Bi,rbpg+α<B~rbpg<Bi,rbp+1g-α; if rbp=1,(ii)
Bi,rbp-1g+α<B~rbpg<Bi,rbp+1g-α; if rbp=2⋯n-1,(iii)
Bi,rbp-1g+α<B~rbpg<Bi,rbpg-α; if rbp=n,whereα is a user-defined parameter that indicates how close the points can be located. To achieve a good fit of the nonlinearity α must be small, relative to the search space on the abscissa axis defined by vmin and vmax. Note that the horizontal boundaries for the endpoints are delimited by their position and their left or right neighbor, respectively. Bounds for mutations on ordinate axis are fixed and equal for all points. This allows each point to move freely throughout the search space on the ordinate axis defined by wmin and wmax. The vertical and horizontal bounds for a selected point are illustrated in Figure 4. When mutation M.2 is required, the selected point can be changed to a new position within the grey rectangle. Details on how to determine vmin, vmaxwmin, and wmax will be given in Section 3.3.Figure 4
Bounds for mutation M.2. Grey area indicatesB~rbpg and B~n+rbpg feasible space.To achieve a fine-tuning of all nonlinearity points, mutations’ aggressiveness can be controlled through standard deviation (17) as in mutation M.1. Note that Δs=wmax-wmin is constant for all mutations over ordinate axis, while for abscissa axis mutations, Δs can be calculated as(20)Δs=Bi,rbp+1g-α-vmin;ifrbp=1Bi,rbp+1g-Bi,rbp-1g-2αifrbp=2⋯n-1vmax-Bi,rbp-1g-αifrbp=n.Mutation M.2 has great potential to explore the searching space. This genetic operation will locate points where there are slope changes. During first generations, it is useful to shape nonlinearity, while in last ones, it allows a refinement. However, when a point is located where there is a slope change, it could be kept in this location until the end of generations, especially when there are abrupt changes in the slope. Because jumps between points are not allowed with mutation M.2, when a point is kept in a place where there is a significant change of slope, one or more points would remain trapped to the left or right of it. This would lead to having redundant points in a segment that would not require so many or worse, to having a segment (curvature) that would not contain enough points. To avoid this drawback, the exploration in the search space is complemented with mutation M.3.Mutation M.3. This genetic operation is designed to concentrate as many points as possible on the curvatures that nonlinearity can have. Therefore, it will be required that each point can be displaced on the abscissa axis by jumping one or more positions of the other points. Let us define a segment as the horizontal space between two consecutive points (so for n points there will be n-1 segments); then a random integer number rs∈1,n-1 will indicate to which segment the selected point will move. The first half of the offspring vector is found using the following expression:(21)B~jg=Bi,rsg+Bi,rs+1g2ifj=rbpBi,jgotherwisewith j=1⋯n. Note that if rbp=rs or rbp=rs-1, the corresponding point will not make a jump but it will be located at the midpoint between its current position and the position of the point on the right or left, respectively. As in the mutation M.2 to prevent points from getting too close together, the α parameter is also used in this mutation; therefore, a jump is conditioned to the space available in the selected segment to accommodate a new point. Minimum space should be 2α. If this condition is not met, rs must be regenerated to randomly search for another segment.To provide a smooth transition between adjacent segments, gene mutation corresponding to the position on the ordinates axis is performed using a quadratic interpolation. To do that, three neighboring points are required. The second half of the offspring vector is found using the following expression:(22)B~jg=B~j-ng2,B~j-ng,1∗k2k1k0ifj=n+rbpBi,jgotherwisewith j=n+1⋯2n. B~j-ng is the v-coordinate of the selected point to mute which can be found with (21); k0, k1, and k2 are the coefficients of the quadratic polynomial Ψ defined by three adjacent points selected once the new v-coordinate of the point that is mutating is known. The three adjacent points can be selected directly when a point has mutated to the first or last segment,(23)Ψ=fBi,1g,Bi,n+1g;Bi,2g,Bi,n+2g;Bi,3g,Bi,n+3g:ifrs=1fBi,ng,Bi,2ng;Bi,n-1g,Bi,2n-1g;Bi,n-2g,Bi,2n-2g:ifrs=n-1,while if the point has mutated to a nonextreme segment, the three adjacent points can be selected using the two points that define that segment plus one on its right or left. For more effective exploration, a random number r3p∈(0,1] is used for selection:(24)Ψ=fBi,rs-1g,Bi,n+rs-1g;Bi,rsg,Bi,n+rsg;Bi,rs+1g,Bi,n+rs+1g:ifr3p≤0.5fBi,rsg,Bi,n+rsg;Bi,rs+1g,Bi,n+rs+1g;Bi,rs+2g,Bi,n+rs+2g:otherwise.Figure5 illustrates how a jump occurs with mutation M.3. Notice that the new ordinate is calculated according to the polynomial formed by the two points of the segment plus a point to the right: that is, r3p>0.5. After a jump has occurred, an ascending reordering of the points with respect to the abscissa values is necessary.Figure 5
Mutation M.3 withrbp=2, rs=4, and r3p>0.5. Jump to the selected segment (dashed line). Quadratic polynomial (solid line).Crossover C.2. This genetic operation works just like crossover C.1 and is applied only to vary the position of a point on the ordinate axis:(25)B~jg=Bi,jg+Bbest,jg2ifj=n+rbpBi,jgotherwisewith j=n+1⋯2n. Bbest,jg is the jth element of vector Bbestg, which corresponds to the individual in the current population g with the best fitness value.


## 3.2.3. Pole-Zero Classification

Due to the stochastic nature of evolutionary algorithms, the binary values of (11) will change as the algorithm evolves, generating different structures of Gw(z) and Gh(z). The evolution of this piece of genetic information is handled by a simple mutation operator.Mutation M.4. Unlike the previous ones, this operator generates a new vector C~g that depends entirely on the effects of mutation, meaning that for this piece of genetic information there is no information exchange between generations. This allows free testing of different structures for Gw(z) and Gh(z) to avoid premature convergence. When this operation is required, a random process will generate the mutation vector:(26)C~jg=1ifNc≤0.50otherwisewith j=1⋯nc+nr+mc+mr. C~jg is the jth element of vector C~g. Nc is a random number with standard uniform distribution on the open interval (0,1). Note that the structure of C~g is built under two considerations: the LTI subsystems cannot be improper and the sum of zeros and the sum of the poles between both subsystems must be equal to the number of zeros and poles of the BLA, respectively.


## 3.3. WH-EA Description

WH-EA is an elitist evolutionary algorithm that evolves a population ofNP individuals. Each individual contains three portions of genetic information related to the parameters of a Wiener-Hammerstein model θig=Pig,Big,Cig. Like any other evolutionary algorithm, WH-EA is inspired by biological evolution over generations. Starting from an initial population, new generations are created using information of the current generation g and performing crossover and/or mutation operations and selection based on the fitness of the new individuals. Algorithm 1 shows a pseudocode of main steps performed in WH-EA, whereas details of all the parameters that the algorithm uses are shown in WH-EA Parameters.Algorithm 1:Pseudocode of WH-EA.
(1) Initialise the population;
(2) Evaluate fitness of all population;
(3) for  g=1  toMaxGendo
(4)Find  θbestg
(5)Random selection of a individual (r1);
(6)Compute γ(g);
(7)if  rpznl≤γ(g)   then
(8)Compute B~g using Algorithm 2;
(9)else
(10)Compute P~g using Algorithm 3;
(11)end if
(12)if  rc≤ξthen
(13)Compute C~g using Mutation M.4;
(14)end if
(15)Update population;
(16) end for
(17) Print θbestMaxGenAlgorithm 2:Modify two-dimensional points for nonlinear function.
(1)if  rnmc≤δnlthen
(2)Compute η(g)
(3)if  rmm≤ηmin+η(g)then
(4)Mutation M.2;
(5)else
(6)Mutation M.3;
(7)end if
(8)else
(9)Crossover C.2;
(10)end ifAlgorithm 3:Modify polo/zero locations.
(1)if  rlmc≤δzpthen
(2)Mutation M.1;
(3)else
(4)Crossover C.1;
(5)end ifInitialize the Population. The initial population P0B0C0 contains NP individuals generated within a search space delimited by lower and upper bounds. For each piece of genetic information, the lower and upper bounds can be determined from the information provided by the BLA (location of poles and zeros and static gain).Real and imaginary values of poles and zeros from the BLA without modifications are introduced as part of the first individualP10 of the initial population. For the rest, mutation M.1 (16) is used but fixing the parent vector as P10, that is, the BLA.Therefore, mutation M.1 must be executedNP-1 times under the above conditions to generate mutated versions of the BLA. Initialization and mutation M.1 are performed considering a search space delimited by lower bound Pmin and upper bound Pmax which are set around the values of the BLA:(27)Pjmin=P1,j0-Υjmin,Pjmax=P1,j0+Υjmax,where Υjmin and Υjmax represent the jth elements of the user-defined vectors Υmin and Υmax, respectively, indicating how much the BLA’s pole/zero position can change as the algorithm evolves.On the other hand, the initial population corresponding to two-dimensional points must be generated within a search space defined horizontally by the minimum (vmin) and maximum (vmax) amplitude of v(t) (the output of Gw(z)) and vertically by the minimum (wmin) and maximum (wmax) amplitude of w(t) (the input to Gh(z)). Although v(t) and w(t) are not known, vmin and vmax depend on the input u(t) and Gw(z). Since Gw(z) is a linear system, the following expressions can be used:(28)vmin=Ω∗umin,vmax=Ω∗umax,where umin and umax are the minimum and maximum values of the signal u(t), respectively, and Ω is a scaling factor which depends on the pole/zero locations and the static gain of Gw(z). Without loss of generality, it is possible to model a Wiener-Hammerstein system considering that both linear blocks have unit gain.Regardingwmin and wmax, if the input u(t) with mean umean enters Gw(z) regardless of its structure, the mean of the output signal vmean will be equal to umean. This is not the case when the signal v(t) enters the nonlinear block since the gain KNL might provide an offset to w(t). However, since Gh(z) is a linear block, the mean of the output ymean will be equal to the mean of the incoming signal wmean. With this information, the straight line of Figure 6 can be drawn, and wmin and wmax can be found using the following equations:(29)wmin=ymean+KNLΩumin-umean,wmax=ymean+KNLΩumax-umean.Figure 6
Information to define the search space for the nonlinear function.As can be seen from (28) to (29), the search space for nonlinearity depends on the input signal, output signal, static gain of the BLA, and Ω which is a user-defined parameter. Since both linear subsystems will be estimated with unit gain, neither of these will amplify their input signals, therefore vmin>umin and vmax<umax. For these two conditions to be met, Ω must be less than one. In the same way, it must be observed that wmin<ymin and wmax>ymax. If Ω is less than one, the first pair of conditions will always be met; however, there is no guarantee that the second pair of conditions will be met. Since it is possible to perform this check prior to the execution of the algorithm, if the second pair of conditions are not met, Ω must be increased, but considering that it must be less than one. It should also be taken into account that if Ω is too large, the search space will be larger than necessary, so the algorithm will cost more to estimate static nonlinearity.To initialize this portion of genetic information, then points are uniformly distributed between Ωumin and Ωumax and located on the straight line shown in Figure 6. These points are introduced as part of the first individual in the population B10. The corresponding genetic information for the rest of individuals is generated using mutation M.2 (19), but considering that the parent vector is always B10.Finally, genetic information corresponding to pole-zero classification is initialized directly using mutation M.4 (26) NP times.Evaluate Fitness. Performance of each individual in the population is defined by a fitness criterion which can be calculated using (14), where θ is obtained from the encoded information in Pig,Big,Cig.The Offspring. Once population has been initialized, for each generation a random integer number r1∈1,NP will be used to select the parent from which an offspring P~g will be generated. As can be seen in Algorithm 1, not all genetic operators are applied at the same time to generate offspring; this can help to expand diversity and avoid premature convergence.One or two pieces of the offspring genetic information will be randomly selected for modification according to their respective genetic operators. A random numberrpznl∈(0,1] chooses between modifying the portion related to static nonlinearity using Algorithm 2 or the portion of genetic information related to pole/zero locations using Algorithm 3. The probability for this selection is handled by the control parameter γ(g) defined as(30)γg=γini1+g∗γrat,γrat=γini/γend2-1MaxGen-1,where γrat is the rate at which the probability γ(g) will decrease from initial probability γini to final probability γend as generations pass; therefore, 0<γend<γini≤1. If γini=1, the probability of modifying the genetic information of nonlinearity in the first generations will be high, while the probability of modifying the location of poles and zeros will be low. On the other hand, if and γend=0.5, in the final generations, the algorithm will modify with equal probability both portions of genetic information. The selection of these values is justified by the fact that pole/zero locations are known and they will only be fine-tuned within a suitable search space to amend possible errors in the BLA estimation, whereas nonlinearity is completely unknown, so the algorithm should focus more on this portion of genetic information during first generations.Variation of genetic information corresponding to the classification of poles and zeros for both LTI subsystems is handled by a comparison between a random numberrc∈(0,1] and the probability ξ∈(0,1]. The value of probability ξ is defined by the user and will be constant throughout the evolution of the algorithm. Figure 7 shows the behaviour of the control parameters (probabilities) used to select the portions of genetic information that will be modified in each generation.Figure 7
Control parameters of WH-EA for selection of the genetic information to be modified in generationg.Algorithm2 is used to modify the genetic information related to nonlinear static function. The control parameter δnl∈(0,1] indicates the probability with which the mutation (either M.2 or M.3) or crossover C.2 will be used. Probability of selecting M.2 or M.3 is variable with respect to the generations. During first generations, mutation M.3 is not necessary, since the nonlinearity can be captured thanks to the two-dimensional points movements due to mutation M.2 and crossover C.2 operations. Since it is very likely that nonlinearity includes one or more curvatures, as the algorithm evolves mutation M.3 will be required to concentrate as many points as possible on these curvatures. The variable probability for selection between both mutations is defined by(31)ηg=1-ηmin-1-ηmingMaxGen,where ηmin∈(0,0.5] is a user-defined parameter indicating the minimum probability with which the mutation M.2 can be selected. Note that according to (31) and Algorithm 2, the maximum probability is 1 and occurs in the first generation. As the algorithm evolves, this probability will decrease linearly until it reaches ηmin in the last generation. When Algorithm 2 is required a random number rnmc∈(0,1] will allow us to select either a mutation or crossover C.2. If a mutation is selected, a new random number rmm∈(0,1] will allow us to select between mutation M.2 or mutation M.3.On the other hand, Algorithm3 is used to modify the genetic information related to pole/zero locations using mutation M.1 or crossover C.1. The control parameter δzp∈(0,1] indicates the probability with which each genetic operation will be used. Since crossover C.1 causes offspring to inherit genetic information from the best individual, a small value of δzp may lead to premature convergence, whereas a value closer to 1 will cause the algorithm to converge very slowly. When Algorithm 3 is required, a random number rlmc∈(0,1] will determine the genetic operation to be used.Update. It is based on a competition between the generated offspring and the individuals of the population. The contestant with the best fitness will be the one who wins the competition. From a randomly selected individual, the offspring starts to compete until defeating an individual; when this happens the descendant will take his place in the population and the algorithm continues with the next generation. If the offspring comes to compete with all individuals and could not win, this will be discarded and the algorithm will pass to the next generation.


## 4. Application of WH-EA and Results

WH-EA was tested on a numerical example and on the benchmark for nonlinear system identification in (SYSID’09) [9], where a Wiener-Hammerstein system is selected as test object. The benchmark is not intended as a competition, but as a tool to compare the possibilities of different methods to deal with this specific nonlinear structure.For both cases, the BLA was estimated with the MATLAB System Identification Toolbox [60] using a Box-Jenkins (BJ) structure. Besides, trends and means were only removed for the BLA identification. The following parameters of the algorithm were set in common for both estimates: ξ=0.25; δzp=0.75; δnl=0.75; ηmin=0.35; in addition, initial and final standard deviations for mutations were set to 20 and 1, respectively.


### 4.1. Numerical Example

A Wiener-Hammerstein system with the following structure was designed (where tansig is the hyperbolic tangent sigmoid transfer function):(32)Gw′z=0.1190z-0.9048,w′t=0.45tansig⁡2.80v′t,Gh′z=-0.01426z-1.0510z+1z-0.9746+0.03656jz-0.9746-0.03656j.A Gaussian excitation signal of 6 dB was filtered with a cut-off frequency of 6 Hz and used as input signal. The system was simulated and 120000 input/output samples were recorded and separated in two parts: the estimation data settn∈1001,70000 for identification purposes and the test data set tn∈71001,120000 for validation purposes (in both data sets first 1000 samples were ignored to avoid transient effects). Furthermore, additive white Gaussian noise with a Signal-to-Noise Ratio (SNR) of 45.32 dB was added to the output.The identification of the BLA was carried out and the model obtained was expressed in factored form:(33)GBLA=-1e-3z-1.0508z+0.9631z-0.9749+0.0366jz-0.9749-0.0366jz-0.9045.The root mean square of the error (eRMS) obtained with this linear model on test data was of 0.0414.According to the BLA structure, vectorP10 was coded with nc=0, nr=2, mc=1, and mr=1 as follows:(34)P10=1.0508,-0.9631,0.9749,0.9045,0.0366.The search space for nonlinear function was defined withΩ=0.28, while elements of Υmin and Υmax were set to 0.01 (except bounds for the zero z=-0.9631 that were set to 0.1, since it influences slightly the dynamics and should have freedom of movement during tuning).WH-EA was executed 3 times withMaxGen=5·106 and different number of points was chosen for the nonlinearity. For all trials, the algorithm was initialized with 60000 individuals and the minimum distance between two points was set to α=(vmax-vmin)/6n. For each estimated Wiener-Hammerstein model, theeRMS on the test data was computed; in addition, the normalized root mean square error (NRMSE) criterion was used to quantify the goodness of fit between real and captured nonlinearity. The results are reported in Table 1; for all cases piecewise linear interpolation was used to connect the n points.Table 1
Performance of the numerical example estimation using different numbers of pointsn to represent static nonlinearity.
n
NRMSE (%)
e
R
M
S
8
99.179
1.183
e
-
3
10
99.396
1.109e-3
12
99.564
1.044e-3The poles and zeros of the BLA were correctly classified in the three tests carried out. As can be seen in Table1, theeRMS of the Wiener-Hammerstein models decreased as the quality of the captured nonlinearity increases. A reasonable model was scored with 12 points considering that the RMS noise was 9.98e-4. A convergence graph for this model is shown in Figure 8, during WH-EA execution, at 6e4th generation the poles and zeros of the BLA were correctly classified, and from there the best individual of each generation conserved the genetic information for this classification. Since the noise RMS is known, at g=1e6 the performance of the model was good enough, so the algorithm could have been stopped. Anyway, 5e6 generations have been allowed in order to demonstrate the great precision that the algorithm can achieve.Figure 8
Convergence graph (a) and NRMSE of the captured nonlinearity (b) for Wiener-Hammerstein estimation withn=12.
(a)
(b)In Figure9, pole/zero locations of the BLA, the obtained Wiener-Hammerstein model, and real system are compared. Notice how WH-EA has moved initial locations trying to get to the true values improving modelling error.Figure 9
(a) Poles of the real system (black +), the BLA (red ×), and the estimatedGw (blue ×) and Gh (green ×) models. (b) Zeros of the real system (black ⋄), the BLA (red o), and the estimated Gh model (green o).
(a)
(b)On the other hand, a graphical comparison between real and captured nonlinearity is shown in Figure10. Linear subsystems of the estimated Wiener-Hammerstein model are represented by (35), while the ordered pairs for the nonlinear static function are shown in Table 2.(35)Gwz=0.0259z-0.9048,Ghz=-0.01960z-1.0512z+0.9713z-0.9746+0.0365jz-0.9746-0.0365j.Table 2
Coordinates of the estimated nonlinearity withn=12.
i
1
2
3
4
5
6
7
8
9
10
11
12
v
i
−0.9694
−0.5794
−0.3716
−0.2479
−0.1251
0.1038
0.1996
0.2893
0.3363
0.4739
0.7648
1.0808
w
i
−0.3303
−0.3250
−0.2894
−0.2354
−0.1415
0.1178
0.2031
0.2558
0.2754
0.3126
0.3309
0.3317Figure 10
Comparison between true (solid-blue) and estimated nonlinearity defined as a piecewise linear function (dashed-red) using pairs[vi,wi] (red circles).


### 4.2. Nonlinear System Identification Benchmark

The system to be modelled is an electronic nonlinear circuit with a Wiener-Hammerstein structure (see Figure11). This system was built by Vandersteen [61] and presented as a benchmark problem for system identification by Schoukens et al. [9].Figure 11
Wiener-Hammerstein benchmark.The first linear dynamic systemG1(s) is designed as a third-order Chebyshev filter (pass-band ripple of 0.5 dB and cut-off frequency of 4.4 kHz). The second linear dynamic system G2(s) is a third-order inverse Chebyshev filter (stop-band attenuation of 40 dB starting at 5 kHz). This system has a transmission zero in the frequency band of interest. This can complicate the identification significantly, because the inversion of such a characteristic is difficult. The system was excited with a filtered Gaussian signal (cut-off frequency 10 kHz). Data used for estimation corresponds to interval tn∈[1,100000], whereas test data corresponds to the remaining part tn∈101001,188000. In order to analyse the performance of estimation methods, the mean value of the simulation error (μ), the standard deviation of the error (std), and the root mean square value of the error (eRMS) must be calculated on test and estimation data [9].Since first 5000 data samples just contain quantization noise, a set of 95000 input/output datatn=5001,…,100000 was used to estimate the BLA. Multiple simulations were performed considering different combinations of poles and zeros for the input/output model and for the noise model. For each BJ model, theeRMS on test data set tn=101001,…,188000 was computed. The BLA was obtained with 6 poles, 5 zeros, and one sample delay for the input/output model and 3 zeros and 3 poles for the noise model. The BLA is fully described with KNL=0.7840 and the pole-zero pattern shown in Figure 12. TheeRMS of this linear model was of 56.159 mV on test data and 43.143 mV after removing trends and means. According to the BLA structure, vector P10 was coded with nc=1, nr=4, mc=2, and mr=2 as follows:(36)P10=0.7605,-0.2733,0,-3.4122,-30.2553,0.6501,0.7314,…,0.8912,0.8289,0.7004,0.4358,0.1692.Figure 12
Poles (×) and zeros (o) of the BLA for the benchmark data (SYSID’09). Two real zeros fall outside the plot in −30.2553 and −3.4122.During BLA estimation stage, different noise models were tested and it was observed that all poles, real zeros within the unitary circle, and complex zeros where they are located correspond to the dominant dynamics of the system, while real zeros outside the unitary circle were more likely to vary their location. This information was used to define the search space for refining the location of poles and zeros:(37)Υmin=0.025,0.025,0.025,1,10,0.025,0.025,…,0.025,0.025,0.025,0.025,0.025,Υmax=0.025,0.025,0.025,1,10,0.025,0.025,…,0.025,0.025,0.025,0.025,0.025.Bounds (37) limit search space in the system dominant dynamics within ±0.025, while for z=-30.255 and z=-3.412, limits are between ±10 and ±1, respectively.Static nonlinearity is represented by piecewise linear functions withn=8 points. Its search space was defined with Ω=0.51 and the minimum distance between two points was calculated with α=(vmax-vmin)/10n. The algorithm was initialized with 5000 individuals and 3e7 generations were executed. The performance of the estimated Wiener-Hammerstein model is shown in Table 3. In Figure 13, it is depicted how the algorithm has distributed pole/zero locations for both linear subsystems. Notice how some of them were displaced to improve the modelling error. The outputs of the estimated Wiener-Hammerstein model, the linear model error, and the nonlinear model error on test data are shown in Figure 14, while the DFT spectra of this signals are shown in Figure 15. Captured nonlinearity is plotted in Figure 16.Table 3
Performance indicators of the estimated Wiener-Hammerstein model. All values are shown in mV.
BLA
Wiener-Hammerstein
Estimation
Test
Estimation
Test
μ
−35.825
−35.951
5.6
e
-
4
1.1
e
-
5
s
t
d
42.108
43.143
0.322
0.306
e
R
M
S
55.286
56.159
0.322
0.306Figure 13
(a) Poles (×) and zeros (o) of the linear subsystemGw (blue) and BLA (red). Zero in z=-30.2553 belonging to BLA and its adjusted final value in z=-26.7609 fall outside the plot. (b) Poles (×) and zeros (o) of the linear subsystem Gh (green) and BLA (red).
(a)
(b)Figure 14
Model output (blue), simulation error of the BLA (grey), and simulation error of the estimated Wiener-Hammerstein model (black).Figure 15
DFT spectra of the modelled output signal (blue), linear model error (grey), and nonlinear model error (black).Figure 16
Captured nonlinearity as a piecewise linear function withn=8 by WH-EA from the benchmark data. Notice that the nonlinear block characterization only needs 14 parameters since the first and last straight segments can be defined just with their angles.Final estimated linear blocksGw(z) and Gh(z) are shown in (38) and (39) respectively, while coordinates for nonlinear function are shown in Table 4. The estimated Wiener-Hammerstein model contains 26 parameters of which 14 are used to represent the static nonlinearity (without end points since they can be located anywhere on their respective end segments; nevertheless, these segments slopes are taken into account). Figure 16 shows how mutation M.2 and mutation M.3 located the two-dimensional points to capture the nonlinearity. As expected, due to the effect of mutation M.3, most of them were concentrated on the curvature.(38)Gwz=6.5e-4z+0.0138z+2.9034z+26.76z-0.7243z-0.7324+0.4361jz-0.7324-0.4361j,(39)Ghz=0.0120z+0.2635z-0.7575+0.6513jz-0.7575-0.6513jz-0.8191z-0.8899+0.1688jz-0.8899-0.1688j.Table 4
Nonlinearity coordinates (n=8, 14 parameters) estimated by WH-EA from benchmark data.
i
1
2
3
4
5
6
7
8
v
i
−0.2168
0.1596
0.3819
0.4943
0.6047
0.7596
1.0811
1.3953
w
i
−0.1979
0.1440
0.3443
0.4276
0.4822
0.5248
0.5698
0.5901


### 4.3. Discussion

In contrast to other methods which generate good initial estimates by splitting the poles and zeros of the BLA, WH-EA allows us to identify Wiener-Hammerstein models avoiding high user interaction which is an advantage compared to methods using QBLA, where at least two intermediate procedures are required before fine-tuning all parameters of the Wiener-Hammerstein model.TheeRMS of 0.306 mV achieved with WH-EA on test data is quite acceptable considering that the RMS of the quantization noise is 0.189 mV. With respect to the initial model (BLA), the error was reduced by a factor of 183.52 thanks to the captured nonlinearity and the updated pole/zero locations. Table 5 shows other proposals that have been tested on the benchmark. It can be appreciated that theeRMS of this paper is slightly higher than others; however, not all estimated models have the same complexity. Some of them use complex models with a greater number of parameters processing raw data before identification, while in this work, WH-EA is fed raw input/output data without preprocessing operations.Table 5
Performance measurements on benchmark data (SYSID’09). All the values are shown in mV.θ indicates the number of parameters used for the model.
Method/technique
e
R
M
S (mV)
θ
Nonparametric BLA, QBLA [41]
0.278
44
Classification of poles and zeros using QBLA [42]
0.286
26
Fractional model parameterization [43]
0.295
26
Advanced method [44, 45]
0.30
64
WH-EA (this paper)
0.306
26
Brute force method [45]
0.31
30
Scanning technique [46]
0.370
-
Polynomial nonlinear state space [47]
0.42
797
Generalized Hammerstein-Wiener [48]
0.481
47
Incremental nonlinear optimization [49]
0.679
25
LS-SVMs [50]
4.070
-
Biosocial culture [51]
8.546
34Comparing WH-EA with the proposals of Westwick and Schoukens [42] and Vanbeylen [43] whose models have the same complexity as the model estimated in this paper, the results are quite similar; however, to obtain a good final model with these two proposals, it is required that the BLA be estimated with high precision. In Vanbeylen [43] at the BLA division phase, a false position of a pole or zero could cause the values of the fractional powers to be close to 1/2 which would cause the user to make a bad decision and the BLA is badly divided. This problem is much more critical in Westwick and Schoukens [42] since the method is based on a graphical comparison between the poles and zeros of the BLA and the poles and zeros of the QBLA. With WH-EA, this problem does not occur, since the evolutionary algorithm contemplates possible errors that can be made in the estimation of the BLA. During the algorithm evolution, the binary code used for the classification of the poles and zeros of the BLA can be changed without user interaction as the false positions of the poles and zeros are corrected. This is an important advantage of WH-EA, since it is very likely that the BLA estimate is subject to errors due to noise and nonlinearity effects; this has been experimentally demonstrated; for this reason, many proposals carry out a final readjustment of the parameters of the Wiener-Hammerstein model.About the computational complexity, an iteration in WH-EA involves the random selection of an individual from the population, the offspring generation according to the genetic operations indicated in Section3, the calculation of the objective function of the offspring, and the population update. This algorithm was implemented on MATLAB® software and was run on a personal computer with core i7 processor of 2,6 GHz and 16 Gb of RAM. The objective function was easily implemented using the filter command to simulate both linear subsystems, while the nonlinearity was interpolated using the interp1 command. For the benchmark system, the average time consumed per 50 iterations was 1.16 s.


## 4.1. Numerical Example

A Wiener-Hammerstein system with the following structure was designed (where tansig is the hyperbolic tangent sigmoid transfer function):(32)Gw′z=0.1190z-0.9048,w′t=0.45tansig⁡2.80v′t,Gh′z=-0.01426z-1.0510z+1z-0.9746+0.03656jz-0.9746-0.03656j.A Gaussian excitation signal of 6 dB was filtered with a cut-off frequency of 6 Hz and used as input signal. The system was simulated and 120000 input/output samples were recorded and separated in two parts: the estimation data settn∈1001,70000 for identification purposes and the test data set tn∈71001,120000 for validation purposes (in both data sets first 1000 samples were ignored to avoid transient effects). Furthermore, additive white Gaussian noise with a Signal-to-Noise Ratio (SNR) of 45.32 dB was added to the output.The identification of the BLA was carried out and the model obtained was expressed in factored form:(33)GBLA=-1e-3z-1.0508z+0.9631z-0.9749+0.0366jz-0.9749-0.0366jz-0.9045.The root mean square of the error (eRMS) obtained with this linear model on test data was of 0.0414.According to the BLA structure, vectorP10 was coded with nc=0, nr=2, mc=1, and mr=1 as follows:(34)P10=1.0508,-0.9631,0.9749,0.9045,0.0366.The search space for nonlinear function was defined withΩ=0.28, while elements of Υmin and Υmax were set to 0.01 (except bounds for the zero z=-0.9631 that were set to 0.1, since it influences slightly the dynamics and should have freedom of movement during tuning).WH-EA was executed 3 times withMaxGen=5·106 and different number of points was chosen for the nonlinearity. For all trials, the algorithm was initialized with 60000 individuals and the minimum distance between two points was set to α=(vmax-vmin)/6n. For each estimated Wiener-Hammerstein model, theeRMS on the test data was computed; in addition, the normalized root mean square error (NRMSE) criterion was used to quantify the goodness of fit between real and captured nonlinearity. The results are reported in Table 1; for all cases piecewise linear interpolation was used to connect the n points.Table 1
Performance of the numerical example estimation using different numbers of pointsn to represent static nonlinearity.
n
NRMSE (%)
e
R
M
S
8
99.179
1.183
e
-
3
10
99.396
1.109e-3
12
99.564
1.044e-3The poles and zeros of the BLA were correctly classified in the three tests carried out. As can be seen in Table1, theeRMS of the Wiener-Hammerstein models decreased as the quality of the captured nonlinearity increases. A reasonable model was scored with 12 points considering that the RMS noise was 9.98e-4. A convergence graph for this model is shown in Figure 8, during WH-EA execution, at 6e4th generation the poles and zeros of the BLA were correctly classified, and from there the best individual of each generation conserved the genetic information for this classification. Since the noise RMS is known, at g=1e6 the performance of the model was good enough, so the algorithm could have been stopped. Anyway, 5e6 generations have been allowed in order to demonstrate the great precision that the algorithm can achieve.Figure 8
Convergence graph (a) and NRMSE of the captured nonlinearity (b) for Wiener-Hammerstein estimation withn=12.
(a)
(b)In Figure9, pole/zero locations of the BLA, the obtained Wiener-Hammerstein model, and real system are compared. Notice how WH-EA has moved initial locations trying to get to the true values improving modelling error.Figure 9
(a) Poles of the real system (black +), the BLA (red ×), and the estimatedGw (blue ×) and Gh (green ×) models. (b) Zeros of the real system (black ⋄), the BLA (red o), and the estimated Gh model (green o).
(a)
(b)On the other hand, a graphical comparison between real and captured nonlinearity is shown in Figure10. Linear subsystems of the estimated Wiener-Hammerstein model are represented by (35), while the ordered pairs for the nonlinear static function are shown in Table 2.(35)Gwz=0.0259z-0.9048,Ghz=-0.01960z-1.0512z+0.9713z-0.9746+0.0365jz-0.9746-0.0365j.Table 2
Coordinates of the estimated nonlinearity withn=12.
i
1
2
3
4
5
6
7
8
9
10
11
12
v
i
−0.9694
−0.5794
−0.3716
−0.2479
−0.1251
0.1038
0.1996
0.2893
0.3363
0.4739
0.7648
1.0808
w
i
−0.3303
−0.3250
−0.2894
−0.2354
−0.1415
0.1178
0.2031
0.2558
0.2754
0.3126
0.3309
0.3317Figure 10
Comparison between true (solid-blue) and estimated nonlinearity defined as a piecewise linear function (dashed-red) using pairs[vi,wi] (red circles).


## 4.2. Nonlinear System Identification Benchmark

The system to be modelled is an electronic nonlinear circuit with a Wiener-Hammerstein structure (see Figure11). This system was built by Vandersteen [61] and presented as a benchmark problem for system identification by Schoukens et al. [9].Figure 11
Wiener-Hammerstein benchmark.The first linear dynamic systemG1(s) is designed as a third-order Chebyshev filter (pass-band ripple of 0.5 dB and cut-off frequency of 4.4 kHz). The second linear dynamic system G2(s) is a third-order inverse Chebyshev filter (stop-band attenuation of 40 dB starting at 5 kHz). This system has a transmission zero in the frequency band of interest. This can complicate the identification significantly, because the inversion of such a characteristic is difficult. The system was excited with a filtered Gaussian signal (cut-off frequency 10 kHz). Data used for estimation corresponds to interval tn∈[1,100000], whereas test data corresponds to the remaining part tn∈101001,188000. In order to analyse the performance of estimation methods, the mean value of the simulation error (μ), the standard deviation of the error (std), and the root mean square value of the error (eRMS) must be calculated on test and estimation data [9].Since first 5000 data samples just contain quantization noise, a set of 95000 input/output datatn=5001,…,100000 was used to estimate the BLA. Multiple simulations were performed considering different combinations of poles and zeros for the input/output model and for the noise model. For each BJ model, theeRMS on test data set tn=101001,…,188000 was computed. The BLA was obtained with 6 poles, 5 zeros, and one sample delay for the input/output model and 3 zeros and 3 poles for the noise model. The BLA is fully described with KNL=0.7840 and the pole-zero pattern shown in Figure 12. TheeRMS of this linear model was of 56.159 mV on test data and 43.143 mV after removing trends and means. According to the BLA structure, vector P10 was coded with nc=1, nr=4, mc=2, and mr=2 as follows:(36)P10=0.7605,-0.2733,0,-3.4122,-30.2553,0.6501,0.7314,…,0.8912,0.8289,0.7004,0.4358,0.1692.Figure 12
Poles (×) and zeros (o) of the BLA for the benchmark data (SYSID’09). Two real zeros fall outside the plot in −30.2553 and −3.4122.During BLA estimation stage, different noise models were tested and it was observed that all poles, real zeros within the unitary circle, and complex zeros where they are located correspond to the dominant dynamics of the system, while real zeros outside the unitary circle were more likely to vary their location. This information was used to define the search space for refining the location of poles and zeros:(37)Υmin=0.025,0.025,0.025,1,10,0.025,0.025,…,0.025,0.025,0.025,0.025,0.025,Υmax=0.025,0.025,0.025,1,10,0.025,0.025,…,0.025,0.025,0.025,0.025,0.025.Bounds (37) limit search space in the system dominant dynamics within ±0.025, while for z=-30.255 and z=-3.412, limits are between ±10 and ±1, respectively.Static nonlinearity is represented by piecewise linear functions withn=8 points. Its search space was defined with Ω=0.51 and the minimum distance between two points was calculated with α=(vmax-vmin)/10n. The algorithm was initialized with 5000 individuals and 3e7 generations were executed. The performance of the estimated Wiener-Hammerstein model is shown in Table 3. In Figure 13, it is depicted how the algorithm has distributed pole/zero locations for both linear subsystems. Notice how some of them were displaced to improve the modelling error. The outputs of the estimated Wiener-Hammerstein model, the linear model error, and the nonlinear model error on test data are shown in Figure 14, while the DFT spectra of this signals are shown in Figure 15. Captured nonlinearity is plotted in Figure 16.Table 3
Performance indicators of the estimated Wiener-Hammerstein model. All values are shown in mV.
BLA
Wiener-Hammerstein
Estimation
Test
Estimation
Test
μ
−35.825
−35.951
5.6
e
-
4
1.1
e
-
5
s
t
d
42.108
43.143
0.322
0.306
e
R
M
S
55.286
56.159
0.322
0.306Figure 13
(a) Poles (×) and zeros (o) of the linear subsystemGw (blue) and BLA (red). Zero in z=-30.2553 belonging to BLA and its adjusted final value in z=-26.7609 fall outside the plot. (b) Poles (×) and zeros (o) of the linear subsystem Gh (green) and BLA (red).
(a)
(b)Figure 14
Model output (blue), simulation error of the BLA (grey), and simulation error of the estimated Wiener-Hammerstein model (black).Figure 15
DFT spectra of the modelled output signal (blue), linear model error (grey), and nonlinear model error (black).Figure 16
Captured nonlinearity as a piecewise linear function withn=8 by WH-EA from the benchmark data. Notice that the nonlinear block characterization only needs 14 parameters since the first and last straight segments can be defined just with their angles.Final estimated linear blocksGw(z) and Gh(z) are shown in (38) and (39) respectively, while coordinates for nonlinear function are shown in Table 4. The estimated Wiener-Hammerstein model contains 26 parameters of which 14 are used to represent the static nonlinearity (without end points since they can be located anywhere on their respective end segments; nevertheless, these segments slopes are taken into account). Figure 16 shows how mutation M.2 and mutation M.3 located the two-dimensional points to capture the nonlinearity. As expected, due to the effect of mutation M.3, most of them were concentrated on the curvature.(38)Gwz=6.5e-4z+0.0138z+2.9034z+26.76z-0.7243z-0.7324+0.4361jz-0.7324-0.4361j,(39)Ghz=0.0120z+0.2635z-0.7575+0.6513jz-0.7575-0.6513jz-0.8191z-0.8899+0.1688jz-0.8899-0.1688j.Table 4
Nonlinearity coordinates (n=8, 14 parameters) estimated by WH-EA from benchmark data.
i
1
2
3
4
5
6
7
8
v
i
−0.2168
0.1596
0.3819
0.4943
0.6047
0.7596
1.0811
1.3953
w
i
−0.1979
0.1440
0.3443
0.4276
0.4822
0.5248
0.5698
0.5901


## 4.3. Discussion

In contrast to other methods which generate good initial estimates by splitting the poles and zeros of the BLA, WH-EA allows us to identify Wiener-Hammerstein models avoiding high user interaction which is an advantage compared to methods using QBLA, where at least two intermediate procedures are required before fine-tuning all parameters of the Wiener-Hammerstein model.TheeRMS of 0.306 mV achieved with WH-EA on test data is quite acceptable considering that the RMS of the quantization noise is 0.189 mV. With respect to the initial model (BLA), the error was reduced by a factor of 183.52 thanks to the captured nonlinearity and the updated pole/zero locations. Table 5 shows other proposals that have been tested on the benchmark. It can be appreciated that theeRMS of this paper is slightly higher than others; however, not all estimated models have the same complexity. Some of them use complex models with a greater number of parameters processing raw data before identification, while in this work, WH-EA is fed raw input/output data without preprocessing operations.Table 5
Performance measurements on benchmark data (SYSID’09). All the values are shown in mV.θ indicates the number of parameters used for the model.
Method/technique
e
R
M
S (mV)
θ
Nonparametric BLA, QBLA [41]
0.278
44
Classification of poles and zeros using QBLA [42]
0.286
26
Fractional model parameterization [43]
0.295
26
Advanced method [44, 45]
0.30
64
WH-EA (this paper)
0.306
26
Brute force method [45]
0.31
30
Scanning technique [46]
0.370
-
Polynomial nonlinear state space [47]
0.42
797
Generalized Hammerstein-Wiener [48]
0.481
47
Incremental nonlinear optimization [49]
0.679
25
LS-SVMs [50]
4.070
-
Biosocial culture [51]
8.546
34Comparing WH-EA with the proposals of Westwick and Schoukens [42] and Vanbeylen [43] whose models have the same complexity as the model estimated in this paper, the results are quite similar; however, to obtain a good final model with these two proposals, it is required that the BLA be estimated with high precision. In Vanbeylen [43] at the BLA division phase, a false position of a pole or zero could cause the values of the fractional powers to be close to 1/2 which would cause the user to make a bad decision and the BLA is badly divided. This problem is much more critical in Westwick and Schoukens [42] since the method is based on a graphical comparison between the poles and zeros of the BLA and the poles and zeros of the QBLA. With WH-EA, this problem does not occur, since the evolutionary algorithm contemplates possible errors that can be made in the estimation of the BLA. During the algorithm evolution, the binary code used for the classification of the poles and zeros of the BLA can be changed without user interaction as the false positions of the poles and zeros are corrected. This is an important advantage of WH-EA, since it is very likely that the BLA estimate is subject to errors due to noise and nonlinearity effects; this has been experimentally demonstrated; for this reason, many proposals carry out a final readjustment of the parameters of the Wiener-Hammerstein model.About the computational complexity, an iteration in WH-EA involves the random selection of an individual from the population, the offspring generation according to the genetic operations indicated in Section3, the calculation of the objective function of the offspring, and the population update. This algorithm was implemented on MATLAB® software and was run on a personal computer with core i7 processor of 2,6 GHz and 16 Gb of RAM. The objective function was easily implemented using the filter command to simulate both linear subsystems, while the nonlinearity was interpolated using the interp1 command. For the benchmark system, the average time consumed per 50 iterations was 1.16 s.


## 5. Conclusions

In this paper, a new method (WH-EA) to identify Wiener-Hammerstein systems in a single step is proposed. The proposal estimates all parameters of the Wiener-Hammerstein model based on a customized evolutionary algorithm (WH-EA). Unlike conventional procedures, WH-EA is able to look for the best BLA split capturing at the same time the process static nonlinearity with high precision, solving a single optimization problem. The algorithm is fed with the estimated BLA and its pole/zero locations are subtly modified within an adequate search space to allow its fine-tuning, while piecewise linear function is used for the nonlinear block. The performance of this approach has been evaluated through a numerical example with a complex static nonlinearity and through the well-known benchmark data (SYSID’09). The results show that it is possible, using WH-EA, to identify a Wiener-Hammerstein system with a good precision in a parametric framework avoiding high user interaction and drawbacks involved in using the QBLA. Further research will be related to WH-EA extension for nonlinear multivariable systems by using a multiobjective optimization approach.



---

*Source: 1753262-2018-02-20.xml*