1. A computerized method of identifying Hammerstein models with known nonlinearity structures using particle swarm optimization, comprising the steps of: (a) generating M input-output data points from a Hammerstein system to be identified, the Hammerstein system having both linear and nonlinear parts, wherein M is a pre-defined integer and the Hammerstein system is defined by y ( k ) = B ( q - 1 ) A ( q - 1 ) x ( k ) + w ( k ) , where k is an integer, the Hammerstein system having an input u(k), y(k) representing the output of the Hammerstein system, w(k) representing measurement noise of the Hammerstein system, and q representing a set of parameters describing the nonlinear part, with q −1 being the unit delay operator, and B(q −1 ) being defined as B(q −1 )=b 0 +b 1 q −1 +. .. +b m q −m and A(q −1 ) being defined as A(q −1 )=1+a 1 q −1 +. .. +a n q −n , where m and n are integers and (m, n) represents the order of the linear part, and x(k) represents a non-measured intermediate variable which is the output of the nonlinear part, x(k) being given by x(k)=f(θ,u(k)), where f(θ,u(k)) is a function representing the nonlinear part and a 1 ,. .. , a n , b 0 ,. .. , b m and θ are unknown parameters to be determined; (b) identifying the Hammerstein system by estimating the unknown parameters a 1 ,. .. , a n , b 0 ,. .. , b m and θ from the input-output data, the estimation comprising the following steps (c) through (h): (c) generating a set of random initial solutions for zeros and poles of the linear part and also for the sets of unknown parameters a 1 ,. .. , a n , b 0 ,. .. , b m and θ of the nonlinear part within a pre-selected range; (d) evaluating a fitness function F for the set of random initial solutions, wherein the fitness function F is given by F = ∑ k = 1 M ( y ( k ) - y ^ ( k ) ) 2 , where y ^ ( k ) = B ^ ( q - 1 ) A ^ ( q - 1 ) x ^ ( k ) , {circumflex over (B)}(q −1 )={circumflex over (b)} 0 +{circumflex over (b)} 1 q −1 +. .. +{circumflex over (b)} m q −m , Â(q −1 )=1+â 1 q −1 +. .. +â n q −n , and {circumflex over (x)}(k)=ĉ 1 u(k)+ĉ 2 u 2 (k)+. .. +ĉ L u L (k), L being an integer defined by representing the non-measured intermediate variable x(k) as x(k)=c 1 u(k)+c 2 u 2 (k)+. .. +c L u L (k), where c 1 , c 2 ,. .. , c L represent a further set of unknown parameters, and â 1 ,. .. , â n , {circumflex over (b)} 0 ,. .. , {circumflex over (b)} m , represent estimates of the unknown parameters a 1 ,. .. , a n , b 0 ,. .. , b m ; (e) minimizing the fitness function F to generate the estimated unknown parameters â 1 ,. .. , â n , {circumflex over (b)} 0 ,. .. , {circumflex over (b)} m and θ, (f) applying particle swarm optimization to the set of solutions by determining the most fit zeros and poles of the linear part calculated in step (d); (g) generating a new set of solutions composed of the most fit zeros and poles of the linear part and with the estimated parameters â 1 ,. .. , â n , {circumflex over (b)} 0 ,. .. , {circumflex over (b)} m and θ generated in step (e); and (h) repeating steps (d) through (g) for a predetermined number of generations.