System identification is defined ausing observed data (input and output) to construct a mathematical model of a plant.There are various algorithms for linear system identification, while non-linear system identification is a growing branch of science. Due to the diversity of non-linear systems structure and behavior, several algorithms have been proposed to model them. Block-oriented models, which are a combination of dynamic linear blocks and non-linear static blocks, are a subset of nonlinear models. Hammerstein--Wiener (H--W) model is a famous subcategory of block oriented models consisting of a dynamic linear block that is serialised with two static non-linear blocks. The importance of H--W model comes from its simplicity and ability to ?t various existing systems and describe their dynamic behaviour . Since 1998, identification of Hermestain-Weiner nonlinear systems have been discussed in lots of litreture. However, most of them, for complexity reduction, were based upon theome simplistic assumptions like: ignor ing the disturbance or measurement noise, modeling of nonlinear blocks with polynomials or other parametric structures, considering the non-linear functions to be continuous, injective or differentiable, the need to apply particular inputs to the system, or the impossibility of generalizing the problem to multi-input-multi-output systems or with color noise. Some strategies also provide non-closed form solutions by appliying numerical methods without a certain variance. Due to the benefit s of using Hemerstein--Wiener model for identification of many real-world systems, first, this problem and literature about it, are discussed in detail in this report. Then, a mathematical structure in the form of the systems of linear equations is defined to fully describe the model, and to reach a generalizable and comprehensive solution. The proposed method do not need to apply special inputs to the system and identified the model with sampled inputs and outputs. The method is designed to consider entered noise and disturbances to a system with multiple inputs and multiple outputs. Considering the greater number of unknowns in identification problem in contrast to number of equations and in order to achieve the unique response with the desired behaviours and to reduce the range of possible solutions, some assumptions are made for solving the problem. These assumptions include the methods used to interpolate non-linear functions, the functions input range partitioning method, or the method of limiting the value or derivative in the nonlinear functions of input and output. One of the used interpolation methods is the Gaussian processes that are applied in combination with the Bayesian method for considering Gaussian disturbances and noise. The proposed method based on constrained optimization has also been used in combination with several nonlinear functions interpolation methods and has been used to identify a typical industrial system. It should be noted that in order to display the efficiency of analytical results, all the methods are presented along with simulation results .