Kajal Kothari, Utkal Mehta,, Vineet Prasad, and Jito Vanualailai

Abstract—The parameter identification of a nonlinear Hammerstein-type process is likely to be complex and challenging due to the existence of significant nonlinearity at the input side. In this paper, a new parameter identification strategy for a block-oriented Hammerstein process is proposed using the Haar wavelet operational matrix (HWOM). To determine all the parameters in the Hammerstein model, a special input excitation is utilized to separate the identification problem of the linear subsystem from the complete nonlinear process. During the first test period, a simple step response data is utilized to estimate the linear subsystem dynamics. Then, the overall system response to sinusoidal input is used to estimate nonlinearity in the process. A single-pole fractional order transfer function with time delay is used to model the linear subsystem. In order to reduce the mathematical complexity resulting from the fractional derivatives of signals, a HWOM based algebraic approach is developed. The proposed method is proven to be simple and robust in the presence of measurement noises. The numerical study illustrates the efficiency of the proposed modeling technique through four different nonlinear processes and results are compared with existing methods.

I. Introduction

IN industrial applications, various advanced processes are inherently nonlinear. The processes with negligible nonlinearity can be treated as linear to reduce complexity; however significant nonlinearity can not be ignored. Several engineering problems are studied widely using two block-oriented models, namely the Hammerstein model and Wiener model.Identification of these models have been an important topic for a long time. The Hammerstein model can represent various nonlinear systems accurately, e.g., fuel cells, chemical processes, magneto-rheological dampers, supercapacitor systems, battery, DC/DC boost converters, aerodynamic models,and biological processes [1], [2]. This highlights the importance of the Hammerstein model and encourages researchers to explore and develop better techniques.

In the Hammerstein model, the input to the linear subsystem is altered by the nonlinearity. Therefore, nonlinearity significantly affects the dynamics of the actual system. It is noteworthy, that nonlinearity elevates computational overhead during identification along with the overall cost, as one has to estimate not only the linear subsystem but also the nonlinearity. Naturally, this leads to the issue of identifying the entire nonlinear model efficiently with minimal computational effort.

In the literature, various techniques have been reported on nonlinear process identification using integer-order (classical)models. The Hammerstein process was estimated using the special test signal in [3] and further extended for Hammerstein-Wiener processes in [4], [5]. The separate block-oriented noniterative relay feedback (SRF) method was illustrated in [6].Mehta and Majhi [7] presented the non-iterative relay feedback(NRF) method to determine the structure prior to the parameters of the Hammerstein model. The Hammerstein system with time delay has been accurately identified using a recursive least squares method in [8]. Recently, a separate block-oriented parameter identification method for Hammerstein systems using least squares was described in [9]. Furthermore, another special input based identification of Hammerstein-Wiener nonlinear system with noise was discussed in [10]. Even though some efficient integer-order techniques have been developed so far, a fractional domain approach is yet to be fully explored for nonlinear process identification.

Initially, Aounet al. [11] introduced fractional calculus based identification techniques for Hammerstein models. The interconnected complex nonlinear processes and their identification problem were addressed in [12]. A fractionalorder model for a thermal system was illustrated in [13] for large variations in temperature. A novel frequency domain approach was developed and presented in [14] whereby Wiener and Hammerstein dynamics were depicted in terms of poles and zeros of the estimated linear approximation, which produced favorable results. However, it was solved with an iterative minimization algorithm which could be caught in local minima. This work was revisited and extended by Giordano and Sj?berg using a time-domain approach to avoid the local minima problem by [15]. Allafiet al. [16]demonstrated identification of the Hammerstein-Wiener model based on a simplified refined instrumental variable method with some prior knowledge. The operational matrix approach using the fractional Taylor basis was developed and utilized for a numerical solution of linear and nonlinear fractional differential equations in [17]. Zhaoet al. [2]presented a parametric identification technique of commensurate fractional-order quasi-linear kinetic battery model (KiBaM). Parameter identification methods of fractional-order chaotic system with time delay was described in [18] whereby an artificial bee colony algorithm was used to solve the multi-dimensional optimization problem. A modified Volterra least mean square (LMS) algorithm for fractional Hammerstein modeling was suggested in [19]. A modified artificial bee colony algorithm was proposed for nonlinear fractional system parameter estimation in [20].Fractional-order multi-input single-output Hammerstein process identification was well illustrated in [21] which utilized the well known genetic algorithm with recursive least squares. Fractional-order stochastic gradient algorithm was developed and discussed for Hammerstein-type ARMAX system identification in [22]. The state-space model approach for fractional Hammerstein process was described in [23].Neural network based nonlinear fractional process modeling verified in real time with wind turbines was depicted in [24].Recently, Wanget al. [1] presented a fully parametric identification (FPI) technique using two algorithms for a fractional Hammerstein process. This FPI method accurately estimates all the parameters of a commensurate model,however, it can not estimate delay. Another technique for Hammerstein process identification using the LMS algorithm and spline interpolation was presented in [25]. A Haar wavelet based parameter estimation technique for fractional nonlinear processes was proposed in [26]. However, this method was not developed for block oriented processes. Furthermore,linear subsystems were considered without time delay.

All the aforementioned fractional approaches do not separate nonlinearity and linear subsystems. The merits of separate identification of linear and nonlinear subsystems make it attractive for researchers. One can have freedom of selecting any one of the available identification techniques for a linear block. Moreover, it is possible to select a non-iterative method for nonlinearity estimation and can improve accuracy due to less number of unknowns to handle at a time. These key points have motivated authors to develop a systematic approach using fractional calculus theory.

In this paper, a parameter identification strategy is presented for a block-oriented Hammerstein-type nonlinear process where the linear part is a continuous-time fractional-order process with time delay. A special test signal constituting unit step and sinusoidal signals, has been applied to the entire nonlinear process for the purpose of identification. Firstly, the unit step response is utilized for the estimation of linear subsystem parameters. It is observed that the step input activates nonlinearity in the form of amplitude changes. In fact, a step response can directly be used for linear subsystem identification without loss of generality. Subsequently, a sinusoidal excitation is utilized to estimate the nonlinear and linear subsystem static gain. Additionally, the HWOM omits the direct calculation of fractional derivatives and transforms complex expressions into a simple algebraic matrix multiplications. This strategy is aimed at precise identification with reduced complexity. The Riemann-Liouville definition(see Appendix) is utilized in this work to handle the fractional derivatives.

The outline of this paper is as follows: Section II introduces the Hammerstein model and Haar wavelet for fractional-order integration. The proposed technique to identify parameters in the linear part and subsequently, nonlinearities is discussed in Section III. Section IV focuses on four numerical examples for performance evaluation and some concluding remarks are provided in Section V.

II. Hammerstain Models and Haar Wavelet

A Hammerstein-type process is represented using blockoriented models, whereby static nonlinearity is succeeded by a dynamic linear subsystem in a cascaded manner as shown in Fig. 1. The inputu(t) and outputy(t) are the measurable quantities whereas the intermediate signalv?(t) is an experimentally immeasurable quantity. It is assumed thatG(s)is a stable process and the static nonlinear functionf(u) is monotonic and crossing at the origin. Let the linear subsystemG(s)be represented by familiar a fractional-order single-pole model as

Fig. 1. Hammerstein-type nonlinear process model.

where coefficients (a0,a1,b0∈R), fractional-order (α ∈R+),and θ represents the input time delay. The purpose of any identification technique is to represent actual system dynamics with more accurate models. The low-order model of any system is highly desirable especially when modeling with a fractional approach. This is due to memory constraints when implementing fractional-orders in real time. Furthermore, it is well-known that many practical systems in the industry are of low-order dynamics, whereby higher-order complicated models are reduced to lower-order compact forms. The advantage of working with a lower-order compact model is the reduction in complexity and computational overhead for identification procedures and the ease in implementing model-based controller design.

Various types of static nonlinearities can be approximated easily using the polynomial form [7]. In this work, for identification purposes, a nonlinear polynomial function is considered up to the fourth order, which is sufficient to represent static memoryless nonlinearities efficiently with reduced complexity. Therefore, the nonlinear block can be illustrated as follows:

However, one can increase the order of the polynomial if higher accuracy is demanded.

A. Concept of Fractional-Order Integration for Delayed Haar Wavelet

In the proposed method, orthogonal Haar wavelet functions are utilized due to high precision, mathematical simplicity,noise immunity and ease of implementation with other standard algorithms. Computationally, Haar wavelets are faster compared to other functions of the wavelet family and can be defined as [27]

whereT fis the total time period,m=2i+k, wherebyi(i≥0)andk( 0 ≤k≤2i) are integers

Now, an arbitrary functionx(t)∈L2[0,T f] can be written in terms of Haar wavelets for the firstMnumber of terms as

Here,XM=[x0,x1,...,xM?1]Tis the Haar coefficient vector andHM(t)=[h0(t),h1(t),...,hM?1(t)]Tis the Haar function vector. If collocation points are considered asti=(2i?1)T f/2M,i=1,2,...,M, theM-square Haar matrix ΦM×Mcan be defined by

The integration of Haar wavelet can be obtained by its multiplication with square matrixThe algebraic expression of fractional-order integration (FOI) of Haar wavelet is written as

whereM-square matrixis called the HWOM of FOI which is computed as given in [27]

whereFαis theM-square generalized operational matrix of FOI for block pulse functions [28] and can be described as

wheref1=1,fq=qα+1?2(q?1)α+1+(q?2)α+1andq=2,3,...,M. Similarly, in the case of delayed input, the considered delayed Haar wavelet functionHM(t?θ) can be characterized as

where θ is the input time delay andZis the Haar wavelet operational matrix of delay which can be written as given in [29]

whereM-square matrixEis the generalized delay operational matrix [28]

Therefore, performing fractional integration of delayed Haar functionHM(t?θ) using (6) and (9) results in

Thus, (IαHM)(t?θ) is thereby simply obtained by matrix multiplication ofHM(t) withandZmatrices.

In the presented work, the fractional calculus and operational matrix based method [29] developed for singleinput single-output process is incorporated and extended to the approximate Hammerstein model.

III. Identification Technique for Hammerstein Models

Complete identification is accomplished using a special test signal in two steps. It is noted that the input and output satisfies the condition ofu(t)∈L2[0,T f] andy(t)∈L2[0,T f].Furthermore, the input signalustepis persistently exciting for the linear subsystem andusineis persistently exciting for overall nonlinear process. First, the unit step input response is utilized for estimation of linear subsystem parameters.Subsequently, the sine wave input response is utilized for parameters of unknown nonlinearity in the second step. The proposed strategy for identification is depicted in Fig. 2.

Fig. 2. Comprehensive nonlinear process identification scheme.

Fig. 3. Utilization of unit step response for linear subsystem without loss of generality.

A. Identification of the Linear Subsystem

Therefore, the substitution ofV(s)=CU(S) in (1) yields

Let,Cb0=B0, therefore

Equation (14) can be simplified as

Expressing (15) in the time-domain yields

Integrating both sides with order αngives

Measured data can now be expressed in terms of the Haar wavelet as

whereYT=[y1,y2,...,yM],UT=[u1,u2,...,uM] and superscriptTdenotes the transpose. Substituting (18) and (19), (17)can be rewritten as

Replacing the fractional integral and time delay by operational matrices yields

Hence

whereIdenotes theM-square identity matrix.

Now, using an optimization approach one can estimate the unknown parameters,a0,a1, α, θ andB0=Cb0. It is to be noted thatb0is multiplied with nonlinear gainC, therefore an independentb0needs to be obtained in the second step with nonlinearity parametersci(i=1,2,3,4).

Remark 1:The model in (1) represents linear fractional first order model which is a very basic model and generally preferred. However, one can use any order model with a greater number of poles and zeros. Consider the generalizednpole,m-zero model as given below.

whereai(i=0,...,n),bj(j=0,...,m) are constants;αi(i=0,...,n), βj(j=0,...,m) are arbitrary real numbers; and θrepresents the time delay. The final output expressiony(t)for the aforementioned transfer function can be written as

B. Identification of the Nonlinearity

The identification of the nonlinear block is accomplished using a sinusoidal response. As mentioned earlier

Ifu(t) andare expressed in terms of Haar wavelets, then

Substituting (28) into (27) results in

The Haar wavelets are piecewise constant and can be expanded in the form of block pulse functions as

where the vector ψM(t) represents block pulse functions. Using (31), the nonlinear terms of (30) can be simplified as given in the literature [26]

Let

Using the definition of block pulse functions, one can write

and substituting the values of (33) and (34), (32) can be further simplified as

which finally yields

Similarly, one can obtain

Now, rewriting all the terms in (30) using (33), (36)–(38)results in

Furthermore, using (24) one can derive the expression of the linear subsystem output for sinusoidal excitation in terms of intermediate simulated input

In the aforementioned expression,b0is unknown, therefore it can be replaced by theb0=B0/C. Note thatB0is a known parameter from the previous step. Therefore,b0andVTcan be replaced using (39), and (40) can now be rewritten as

This expression of outputy(t) is used to findci(i=1,2,3,4),keeping previously estimated linear subsystem parameters constant. An independentb0is obtained usingb0=B0/Cfrom estimatedcivalues.

The proposed approach is summarized into the following simple steps.

Step 1:Generate a special input signal, constituting a unit step and sinusoidal signal, for excitation of the Hammerstein process.

Step 2:Record the measurable input and output datau(t)andy(t) for 2Msamples.

Step 3:Identify the linear subsystem parametersa0,a1, α, θ and alsoB0fromMsamples (here, a half-time test data is basically recorded for a step input).

Step 4:Calculateb0andcifrom the remaining half-time (M)samples.

The time-moment weighted integral performance criteria,for identification of both linear and nonlinear subsystems, is used as an objective function to minimize identification error.The integral of squared-time-weighted-error (ISTE) is used to estimate the system parameters and can be written in accordance to estimated values

where ρ is the vector of unknown parameters,y(k) andydata(k)are the simulated response and collected actual data respectively, andMdenotes the total number of samples. The objective of the optimization is to find model parameters that would ideally reduce theISTEto minimum value. The MATLAB function fsolve is used to compute the best-estimated parameters which would satisfy the objective function in (42). The function solve has three different algorithms: trust-regiondogleg, trust-region, and Levenberg-Marquardt. In this work,the Levenberg-Marquardt algorithm was utilized and its detailed convergence analysis is depicted in [30]. Moreover, the persistence of excitation is a sufficient condition for parameter convergence [31], and as mentioned earlier, inputs for the proposed method are persistently exciting. Furthermore, owing to the convex nature and feasible solution of problem (42),the convergence of the prevailing mathematical optimization algorithms is always guaranteed. However, it is always important to assume the appropriate initial value of parameters.More appropriate initial guesses make convergence occur more quickly with more accuracy. There is no mathematical formula or theory available for the selection of initial guesses.The proposed identification method gave satisfactory results in the finite trials. The first optimization procedure was performed, and the result was used as the new initial guess after which the succeeding trial is performed and so on [28]. The selection of initial guesses affects the convergence and number of trials required for the optimization procedure.

IV. Numerical Study

The proposed identification scheme for fractional-order Hammerstein models, was implemented and verified in MATLAB. Several numerical examples, studied in the literature in [1], [6], [7], and [11], were simulated to evaluate the performance with a significant nonlinearity. The optimum number of samples is considered on the basis of the acceptable accuracy of identification. The performance was verified in the following examples with respect to number of samples,different input types, and measurement noises.

Example 1:Consider a fractional-order Hammerstein processp1studied in [11] whereby the linear subsystem isG1(s)=e?s/1+2s0.5and nonlinearity isf1(u)=u(t)+2u2(t).The importance of nonlinear process identification becomes apparent when analyzing the function characteristic of a nonlinear element with the 3D plot as shown in Fig. 4. The input compounded of step and sinusoidal signals for the process is illustrated in Fig. 5.

Fig. 4. Physical interpretation of nonlinearity via 3D plot for f1(u).

Fig. 5. Input output data for identification.

As outlined in the previous section, a unit step responseu(t)=Cv?(t)was considered without loss of generality and after that, the linear subsystem parameters,B0,a0,a1, α andθ were estimated from the measured output. Using the second part of measured data from the sinusoidal response, the nonlinearity polynomial coefficients,c1,c2,c3,c4, and independent linear subsystem gainb0, were obtained. Table I shows the approximated linear subsystems, nonlinearity and time-domain errors for studied examples. Figs. 6 and 7 compare the actual and estimated linear and nonlinear parameters.

Noise is a key factor in determining the effectiveness of the technique for any real process identification. The identification algorithm should be robust in order to process real-time measurement values. The proposed method was verified under measurement noises with SNR = 10 dB and SNR = 20 dB. The results plotted in Fig. 8 depict robustness when the measured data is contaminated by noise. This result can be credited to the noise immunity property of the Haar wavelets.

Fig. 6. Frequency response Nyquist plot for G 1(jω).

Fig. 7. Identified nonlinearity for Example 1.

Fig. 8. Results with noisy output signals for Example 1.

The data length is another important factor to determine the trade-off between the accuracy and speed of the identification.Furthermore, the complexity in separately identifying the error introduced is profound due to the operational matrix approximation. Nonetheless, the error due to approximation is dependent on the data lengthM. It is obvious that a small number of samples accelerates the identification but can result in incorrect estimation. Table II shows the various datalengths with time-domain errors. It clearly indicates the effect of the data length on the precision of the identification.Overall, it is observed that the Haar wavelet based method works favorably with sample pointsM=256.

TABLE I Results and Comparisons

TABLE II Effect of Datalength on Accuracy of Identification

In order to verify identification results with other deterministic or random signals, the same processp1is identified with various input types. It is noteworthy that the estimation procedure was considered with the same initial guess, samples, and number of iterations. The resulting timedomain errors are listed in Table III. It is seen that the presented technique is almost consistent with any input test signal. Furthermore, the normalized mean values of all parameters are depicted in Fig. 9 for various inputs.

TABLE III Various Inputs for b0 and Nonlinearity Identification

Fig. 9. Effect of input types on parameter identification for G 1(s).

Example 3:The nonlinear processp3was considered from[6], with linear subsystemG3(s)=e?s/(s+1)5and input nonlinearityf3(u)=(1?e?0.75u)|u|. Using the proposed method, the identification results are shown in Table I. The identification errors and the comparative plots in Figs. 13 and 14 clearly demonstrate the superiority of the proposed method over the SRF [6] method.

Example 4:The nonlinear processp4was considered from the recent literature [1], with linear subsystemG4(s)=1.2/(s3+3s1.5+2) and input nonlinearityf4(u)=u+0.5u2+0.25u3. In order to make a fair comparison, the linear and nonlinear subsystems were estimated with the same number of unknown parameters as in [1]. The proposed method has been enhanced as given in Remark 1 for comparison purposes. The identified linear and nonlinear subsystems with measurement noise SNR = 20 dB are shown in Table I. The identification errors demonstrate similar accuracy. The main advantage of the proposed method over the FPI method [1] is that the proposed method can accurately estimate the system with time delay and non-commensurate orders without prior knowledge.

Fig. 10. Comparison of step responses for linear subsystem G 2(s).

Fig. 11. Comparison of static nonlinearities f2(u).

The convergence analysis for linear and nonlinear subsystems are depicted in Figs. 15 and 16, respectively. The convergence of objective functionISTEcan be achieved towards a minimum possible value in a finite number of trials.The first optimization procedure was performed and the result was used as the new initial guess after which the succeeding trial was performed and so on. It is obvious that for each trial,ISTEvalue decreases as iteration increases. Finally, the optimum values of parameters were obtained with minimum ISTE as shown in Fig. 15.

Fig. 12. Comparison of sinusoidal responses for process p2.

Fig. 13. Comparison of sinusoidal responses for process p3.

Fig. 14. Comparison of static nonlinearities f3(u).

Fig. 15. Convergence analysis of linear subsystem using step response.

Fig. 16. Convergence analysis of nonlinear subsystem using sine response.

V. Conclusions

As seen from the results presented, the Haar wavelet-based identification strategy for a class of the Hammerstein model definitely reduces the complexity with fractional-order estimation. No prior information is needed and it simultaneously identifies the delay parameter with other linear subsystem parameters. It is noteworthy that the presented approach completely separates the linear and nonlinear subsystem identifications. The proposed technique permits the freedom to select any of the available fractional identification techniques for linear models. Moreover, it is possible to select a non-iterative method for a nonlinearity estimation. It is observed from the results that the scheme is robust with noisy signal data, eliminating additional denoising steps. Both linear and nonlinear subsystems are accurately identified using a simple HWOM based algebraic approach. Furthermore, the approach can be enhanced with a higher-order fractional model in order to achieve more accuracy but requires more computational efforts. It would be interesting to further extend this work for the Wiener or Hammerstein-Wiener models.

Appendix Fractional Calculus

Fractional calculus is a generalization of non-integer (real)order integration and differentiation and its operator is generally defined as

whereaandtare the bounds of the operation and α (α ∈R) is the order of operation.

There exist numerous definitions to characterize fractional integration and differentiation. In the proposed work, the Riemann-Liouville definition is utilized which can be written as

wheren?1<α

wheresαis a fractional Laplacian operator. The fractional R-L integration of an arbitrary functionf(t) is given by

By using the convolution property, (46) can be further simplified as

wheret>aand ? denotes convolution. It can be written in the Laplace domain with zero initial conditions as