Instrument variable method based on nonlinear transformed instruments for Hammerstein system identification

Abstract

An extended instrumental variable (EIV) method is considered for the stochastic Hammerstein system (ARMAX and general model structure). The EIV method provides consistent parameter estimates by eliminating noise-induced bias in the least square (LS) method. To estimate the parameters, the Hammerstein model is formulated using the bilinear parameterization. The bilinear model is identified by introducing the nonlinear instrumental variables obtained from transformed delayed outputs using nonlinear mapping and polynomial basis of delayed inputs. These instruments are analyzed in full generality by computing the bounds on expected relationship between instruments and noise for the general noise disturbance structure. Then, a specific case with hyperbolic tangent (tanh) transformation is considered. Comparative performance analysis of the proposed IV method with the existing IV method, the data filtering-based LS methods, and the extended LS method shows improvement in the statistical properties of parameters estimates.

Keywords

Parameter estimation nonlinear instrumental variables least-square method instrumental variable method Hammerstein models stochastic models

1. Introduction

Identification of nonlinear systems has fundamental importance, since a large class of nonlinear processes exhibit nonlinearities that cannot be linearized. The Hammerstein model structure can effectively represent such processes and is composed of a static nonlinearity followed by a linear dynamical system. These models have been successfully used in modeling of several industrial applications, e.g. electrical drives (Balestrino et al., 2001), sticky control valves (Srinivasan et al., 2005), solid oxide fuel cells (Jurado, 2006), stretch reflex dynamics (Westwick and Kearney, 2001), and control design (Kuo et al., 2012).

Several approaches have been studied for the identification of Hammerstein models, such as maximum likelihood estimation (Li et al., 2014), separable least square methods (Bruls et al., 1999; Westwick and Kearney, 2001), stochastic gradient (SG) methods (Xie and Yang, 2010), and the least square (LS) methods (Bai and Chan, 2008; Wang and Ji, 2014). The LS method due to its simple structure has been successfully used for the discontinuous and piecewise nonlinearities for Hammerstein ARX models (Vörös, 1999; Hu and Ding, 2013). The LS estimates despite being convergent are inconsistent due to an asymptotic bias for Hammerstein OE and ARMAX models. To cope with the bias problem, Ding and Chen (2005) proposed the extended LS (ELS) method for Hammerstein ARMAX models, Wang and Ding (2016a) studied the extended SG algorithm for multiple-input multiple-output (MIMO) Hammerstein CARMA models, and Ding et al. (2007) presented the auxiliary model based LS method for OE models. The hierarchical identification methods have been studied for Hammerstein dual-rate models using the polynomial transformation (Wang et. al, 2015) and MIMO Hammerstein systems (Wang, 2016). Recently, a data filtering-based LS (DFLS) algorithm was proposed for multivariable linear time invariant (LTI) systems to improve the estimation efficiency (Wang and Ding, 2016b) and also extended for Hammerstein and Hammerstein–Wiener systems with colored noise (Wang et. al., 2013; Wang and Ding, 2015, 2016c).

Another approach for unbiased estimation is to consider the instrumental variable (IV) method. This approach has been commonly used for discrete LTI systems (Young, 1970; Soderstrom and Stoica, 1981, 1983; Stoica and Soderstrom, 1983). However, certain restrictions need to be imposed on the noise model, e.g. instruments using delayed outputs give consistent estimates when the combined plant and noise belong to the ARMAX model structure (Soderstrom and Stoica, 1981). For Hammerstein Box–Jenkins models, Laurain et al. (2008) proposed the refined IV method in which the instruments are selected by time invariant filtering of observed variables and augmented model outputs. In particular, the selection of appropriate instruments relies on the linear transformation of observed system variables. These observations motivate the authors to explore the IV method using the instruments chosen as delayed outputs and provide an IV method using the nonlinear transformation of delayed outputs to improve the statistical properties of estimates.

The objective of this paper is to study the IV identification of the Hammerstein model, with an emphasis on addressing two problems. The first is to construct the instruments by means of observed delayed outputs and polynomial basis of inputs, and then extend the case by providing the instruments from the nonlinear transformation of delayed outputs. The second problem is concerned with the consistency of the transformed instruments to cover the model structure other than ARMAX. It is shown that the proposed method gives consistent estimates even in the case of noise modeling error with improvement in the statistical properties. Furthermore, this method optimizes the computational load, as it utilizes the instruments which have neither the operation of prefiltering nor the augmented model (by inclusion of noise dynamics). These significant properties make the presented IV method as an efficient method to estimate the parameters of Hammerstein ARMAX and other model structures.

The paper is organized as follows. Section 2 describes the problem formulation. Section 3 discusses the extended IV (EIV) method. The nonlinear transformed instruments are discussed in full generality along with certain conditions in Section 4. Section 5 deals with the consistency and normality conditions. In Section 6, it is implemented and compared with the LS method, the nonlinear LS method, the ELS method, the DFLS methods and the existing IV method through relevant Monte Carlo simulations. Finally, conclusions are drawn in Section 7.

In this paper, $ℝ$ and $ℤ$ ⁺ denotes the set of all real numbers and positive integers respectively. $ℝ$ ⁿ denotes the vector space of all the n-tuples over $ℝ$ (its elements are called n-dimensional vectors). $ℝ$ ⁿ ^× ^m denotes the set of all $n \times m$ matrices over $ℝ$ . Finally, I_m denotes the Identity matrices of order m. The superscript $T$ denotes the matrix transpose. For the symmetric matrix $M \in ℝ n \times n, M ≻ 0$ means that M is positive definite.

2. Hammerstein system model

Consider a single input–single output, Hammerstein system (Figure 1) described as

y (k) = \frac{B (q - 1)}{A (q - 1)} f (u (k)) + η (k)

(1)

η (k) = \frac{C (q - 1)}{A (q - 1) D (q - 1)} ɛ (k)

(2)

where,

u (k)

y (k)

, and

η (k)

are the system input, output, and correlated process generated from white noise

ɛ (k)

, and

A (q - 1), B (q - 1)

C (q - 1), and D (q - 1)

are polynomials defined as

A (q - 1) = 1 + a_{1} q - 1 + a_{2} q - 2 + \dots + a_{n_{a}} q - n_{a} B (q - 1) = b_{1} q - 1 + b_{2} q - 2 + \dots + b_{n_{b}} q - n_{b} C (q - 1) = 1 + c_{1} q - 1 + c_{2} q - 2 + \dots + c_{n_{c}} q - n_{c} D (q - 1) = 1 + d_{1} q - 1 + d_{2} q - 2 + \dots + d_{n_{d}} q - n_{d}

(3)

By setting $D (q - 1) = 1$ in equation (2), the Hammerstein ARMAX structure can be obtained. In order to proceed, $y (k)$ needs to be defined in recursive form as follows

y (k) = \sum_{n = 1}^{n_{a}} - a_{n} y (k - n) + \sum_{m = 1}^{n_{b}} b_{m} f (u (k - m)) + \sum_{j = 0}^{n_{c}} c_{j} ɛ (k - j) : (c_{0} = 1)

(4)

The static nonlinearity $f (u (k)) = \sum {no}_{l = 1}^{n_{p}} p_{l} u l (k),$ is approximated using basis polynomial functions with a finite number of parameters (Stone, 1948). Substituting $f (u (k))$ in equation (4), we get

y (k) = z T (k) θ + e (k)

(5)

e (k) = \sum_{j = 0}^{n_{c}} c_{j} ɛ (k - j) : (c_{0} = 1)

(6)

where

and

θ = [a_{1}, \dots a_{n_{a}}, b_{1} p_{1}, b_{1} p_{2}, \dots b_{1} p_{n_{p}}, \dots, b_{n_{b}} p_{1}, \times b_{n_{b}} p_{2}, \dots b_{n_{b}} p_{n_{p}}] T \in ℝ (n_{a} + n_{b} n_{p})

(7b)

The orders n_a, n_b, and n_p are assumed to be known a priori. An examination of equation (5) reveals that the resulting bilinear model is not unique in terms of parameters. In order to get a unique parameterization, the coefficient of first basis polynomial function is assumed to be 1, i.e. $p_{1} = 1$ (Narendra and Gallman, 1966; Laurain et al., 2008).

For N measurements of input–output data pairs (u(k),y(k)) for $k \in S \subseteq ℤ +, Card (S) = N$ .

where $N >> n_{a} + n_{b} n_{p},$ equation (5) yields the following matrix equation

Y = Z θ + E

(8)

where

Y = [y (1), y (2), \dots y (N)] T \in ℝ N

E = [e (1), e (2), \dots e (N)] T \in ℝ N

and

Z

is a

N \times (n_{a} + n_{b} n_{p})

matrix, whose

k th

row is the vector

z T (k)

To ensure the consistency of instruments, the following standard assumptions are made:

$A_{1}$ . The signals $u (k) and ɛ (k)$ are stationary, ergodic, and have zero-mean so that the ensemble average may be replaced by time averages over one sample function.

$A_{2}$ . The signals $u (k)$ and $ɛ (k)$ are assumed to be statistically independent of each other.

3. The EIV method

Let a vector $q (k) \in ℝ n_{q}$ with dimension $n_{q} \geq n_{a} + n_{b} n_{p}$ (called the vector of instruments) be given such that the following properties hold ( $E []$ – expectation operator)

{lim}_{N \to \infty} \frac{1}{N} \sum_{k = 1}^{N} q (k) e (k) = {lim}_{N \to \infty} {\hat{r}}_{qe} (N) Δ = r_{qe} = E [q (k) e (k)] = 0, w . p . 1

(9a)

{lim}_{N \to \infty} \frac{1}{N} \sum_{k = 1}^{N} q (k) ρ (k) = {lim}_{N \to \infty} {\hat{r}}_{q ρ} (N) Δ = r_{q ρ} = E [q (k) ρ (k)] = 0, w . p . 1

(9b)

where the residual vector

ρ (k)

is defined as

ρ (k) : = y (k) - z T (k) \hat{θ} (N)

. Let the matrix

{\hat{R}}_{qz} = {\hat{R}}_{qz} (N)

and the vector

{\hat{r}}_{qy} = {\hat{r}}_{qy} (N)

be defined in asymptotic terms as

\frac{{lim}_{N \to \infty} \frac{1}{N} \sum_{k = 1}^{N} q (k) z T (k) = {lim}_{N \to \infty} {\hat{R}}_{qz} (N) Δ}{= R_{qz} = E [q (k) z T (k)], w . p . 1}

(10a)

\frac{{lim}_{N \to \infty} \frac{1}{N} \sum_{k = 1}^{N} q (k) y (k) = {lim}_{N \to \infty} {\hat{r}}_{qy} (N) Δ}{= r_{qy} = E [q (k) y (k)], w . p . 1}

(10b)

Now if for a square matrix $W ≻ 0$ , the condition of equation (9) is satisfied and

rank (R_{qz}^{T} {WR}_{qz}) = n_{a} + n_{b} n_{p}, w . p . 1

(11)

then an EIV estimator of

θ

is given by (Stoica and Soderstrum, 1983)

\hat{θ} IV = ({\hat{R}}_{qz}^{T} W {\hat{R}}_{qz}) - 1 {\hat{R}}_{qz}^{T} W {\hat{r}}_{qy}

(12)

Here, $\hat{θ} IV$ represents the IV estimates whereas the true parameters will be denoted by $θ \circ$ . The solution of equation (12) is a constrained solution with constraints $W ≻ 0$ and $p_{1} = 1$ . For convenience, we consider the case $W = I_{n_{q}}$ as it will not affect the consistency of estimates.

Remark 1

The parameter estimates of the Hammerstein model are calculated using the averaging principal. Let $\hat{θ} IV (i)$ be the ith component of $\hat{θ} IV$ , then referring to equations (7b) and (12), we have

{\hat{p}}_{j} = \frac{1}{n_{b}} \sum_{i = 1}^{n_{b}} \frac{\hat{θ} IV (n_{a} + (i - 1) n_{p} + j)}{{\hat{b}}_{i}}, j = 2, 3, \dots, n_{p}

(13)

where

{\hat{b}}_{i} = \hat{θ} IV (n_{a} + (i - 1) n_{p} + 1), 1 \leq i \leq n_{b}

. The estimates and true parameters are represented by

{\hat{θ}}_{s}^{IV} = [{\hat{a}}_{1}, \dots {\hat{a}}_{n_{a}}, {\hat{b}}_{1}, \dots {\hat{b}}_{n_{b}}, {\hat{p}}_{2}, \dots {\hat{p}}_{n_{p}}] T \in ℝ n_{s}

(n_{s} = n_{a} + n_{b} + n_{p} - 1)

and

θ_{s} \circ = [a_{1}, \dots a_{n_{a}}, b_{1}, \dots b_{n_{b}}, p_{2}, \dots p_{n_{p}}] T \in ℝ n_{s},

respectively.

4. The IV method using nonlinear instruments for the Hammerstein model

In this section, instruments obtained using the nonlinear transformation of delayed observed variables for the identification of Hammerstein ARMAX and general model structure are analyzed. In the subsequent subsection, a specific case using the hyperbolic tangent transformation of delayed outputs is discussed along with the adequacy of the fit.

4.1. Nonlinear instruments

The asymptotic efficiency of the IV method depends on the instruments therefore the main goal is to choose appropriate instruments. In this paper, we are interested in the nonlinear transformation $T : ℝ \to ℝ$ of the observed variables (specifically observed output) to form the IV variables which satisfies the following.

$T$ is continuous over $ℝ$ and is bounded in all bounded intervals in $ℝ$ .

$T (x) = O (e | x | α)$ as $| x | \to \infty$ where $α \in ℝ$ and $α < 2$ .

Let $T (y (k) n_{qy}) = g_{i} (y (k), n_{qy})$ , where $g_{i} (y (k), n_{qy}) = g (y (k - i - n_{qy}))$ $, i = 1, \dots,, n_{a}^{'}$ with $n_{a}^{'} \geq n_{a}$ ( $n_{qy}$ represents the delay of observed outputs and $n_{a}^{'}$ represents the number of functions g_i). A large delay may lead to the ill conditioned matrix ${\hat{R}}_{qz}^{T} {\hat{R}}_{qz}$ , resulting in poor estimation accuracy (Soderstrom and Stoica, 1981). The remedy is to use the extended vector such that the condition ${\hat{R}}_{qz}^{T} {\hat{R}}_{qz} ≻ 0$ can be fulfilled easily or at least have better conditioning. Then the IV vector of dimension $n_{q} \geq n_{a} + n_{b} n_{p}$ is given by

q T (k) = [g_{1} (y (k), n_{qy}) \dots, g_{n_{a}^{'}} (y (k), n_{qy}) u (k - 1), u 2 (k - 1), \dots u n_{p} (k - 1), \dots u (k - n_{b}), u 2 (k - n_{b}), \dots u n_{p} (k - n_{b})]

(14)

The method will be abbreviated as Hammerstein nonlinear extended IV method (HNEIV). Suppose that the polynomial basis functions (see equation (5)) are bounded in the given experimental range. This condition is satisfied when $u (k)$ is bounded stationary random sequence. The assumption A₂ together with equation (5) implies that $u (k)$ and $e (k)$ are independent and therefore their functions, namely ϕ( $u (k)$ ) and ψ( $e (k)$ ), are also independent (Grimmett and Stirzaker, 2001). In particular, if ψ is the identity function, then $ϕ (u (k)) = u l (k - m)$ and $ψ (e (k)) = e (k)$ are independent for $l = 1, 2, \dots, n_{p}$ , $m = 1, 2, \dots, n_{b}$ , $k > m$ .

Lemma 1

Suppose the assumptions $A_{1}$ and $A_{2}$ hold. Then

E [g_{i} (y (k), n_{qy}) e (k)] = 0, n_{qy} = n_{c} + 1, n_{c} + 2, \dots . and i = 1, 2, \dots n_{a}^{'}

Proof

Let $r_{y ɛ} (- t) = E [y (k) ɛ (k + t)]$ for $t \in ℤ +$ . Then by the assumption on $ɛ (k)$ , $r_{y ɛ} (- t) = 0$ for all $t \in ℤ +$ . Let $r_{ɛ e} (- t) = E [ɛ (k) e (k + t)]$ . Now by using equation (6), we get

r_{ɛ e} (- t) = E [ɛ (k) ɛ (k + t)] + c_{1} E [ɛ (k) ɛ (k + t - 1)] + \dots + c_{n_{c}} E [ɛ (k) ɛ (k + t - n_{c})]

(15)

Again the assumption that

ɛ (k)

is a white noise sequence immediately leads to

r_{ɛ e} (- t) \neq 0,

t \leq

n_c. Thus

r_{ɛ e} (- t) = 0

t = n_{c} + 1, n_{c} + 2, \dots

Similarly,

r_{ee} (- t) = 0

t = n_{c} + 1, n_{c} + 2, \dots

Now using (4) and

f (u (k)) = \sum {no}_{l = 1}^{n_{p}} p_{l} u l (k)

, we have

r_{ye} (- t) = E [(\sum_{n = 1}^{n_{a}} - a_{n} y (k - n) + \sum_{m = 1}^{n_{b}} b_{m} \times \sum_{l = 1}^{n_{p}} p_{l} (u l (k - m)) + e (k)) e (k + t)]

(16)

By assumption A₂, the functions $u l (k - m)$ where $m = 1, 2, \dots, n_{b}$ and $l = 1, 2, \dots, n_{p}$ are independent of $e (k + t)$ for $t = n_{c} + 1, n_{c} + 2, \dots$ . Consequently $y (k)$ is independent of $e (k + t)$ , and therefore $r_{ye} (- t) = 0$ . Now it follows easily that

E [g_{i} (y (k), n_{qy}) e (k)] = 0 forall n_{qy} = n_{c} + 1, n_{c} + 2, \dots . and i = 1, 2, \dots n_{a}^{'}

▪

4.2. Noise modeling error

In some instances, the ARMAX model cannot describe the power spectrum of noise sufficiently. Therefore, in order to avoid the possible bias, the noise model structure may be modified with $D (q - 1) \neq 1$ (for details see Soderstrom and Stoica, 1989; Ljung, 1999; Schoukens et al., 2007). Then, the noise sequence $e (k)$ can be written as

e (k) = C (q - 1) v (k), v (k) = \frac{1}{D (q - 1)} ɛ (k)

(17)

where

v (k)

is an autoregressive (AR

(n_{d})

) process. Assume that

D (z)

(z being an arbitrary complex variable replacing

q - 1

) has all its zeroes strictly outside the unit circle. Then,

v (k)

can be expressed as infinite parameters moving average (MA(∞)) process (Berk, 1974)

v (k) = F (q - 1) ɛ (k) = \sum_{j = 0}^{\infty} f_{j} ɛ (k - j), (f_{0} = 1), \sum_{j = 0}^{\infty} j β | f_{j} | < \infty

(18)

where

β \geq 0

. Let

\tilde{v} (k)

be the fitted AR

(n_{\tilde{d}})

process based on N observations

v (1), \dots, v (N)

. The process

\tilde{v} (k)

can be written as an invertible MA(∞) process

\tilde{v} (k) = \tilde{F} (q - 1) ɛ (k) = \sum_{j = 0}^{\infty} {\tilde{f}}_{j} ɛ (k - j), ({\tilde{f}}_{0} = 1)

(19)

such that

sup {no}_{N > N_{1}} \sum {no}_{j = 0}^{\infty} j β | {\tilde{f}}_{j} | < \infty

for some

N_{1} \in ℤ +

(Buhlmann, 1995). The coefficients f_j and

{\tilde{f}}_{j}

\forall j \in ℤ +

satisfy

sup {no}_{j \in ℤ +} | {\tilde{f}}_{j} - f_{j} | = O (\sqrt{(logN / N)})

. Due to the stationarity,

{\tilde{f}}_{j} \forall j \in ℤ +

dies out exponentially as j becomes large (Box and Pierce, 1970). This dying off is so fast that the summation can generally be stopped at a value of j much less than N, unless

D (z)

is extremely close to the boundary of stationary region. Here we assume that there exists

n_{\hat{f}} \in ℤ +

such that

sup {no}_{j > n_{\hat{f}}} | {\tilde{f}}_{j} | = O (1 / \sqrt{N})

. Then, the truncated MA

(n_{\hat{f}})

process is given by

\hat{v} (k) = \hat{F} (q - 1) ɛ (k) = \sum_{j = 0}^{n_{\hat{f}}} {\hat{f}}_{j} ɛ (k - j), ({\hat{f}}_{0} = 1) where {\hat{f}}_{j} = {\tilde{f}}_{j} forall j \leq n_{\hat{f}}

(20)

Lemma 2

Let $F (z) = \sum {no}_{j = 0}^{\infty} f_{j} z j$ (z being an arbitrary complex variable replacing $q - 1$ ) be bounded away from zero for $| z | \leq 1$ and $\sum {no}_{j = 0}^{\infty} j r | f_{j} | < \infty$ . Suppose $sup {no}_{j \in ℤ +} | {\tilde{f}}_{j} - f_{j} | = O (\sqrt{logN / N}) and$ $sup {no}_{j > n_{\hat{f}}} | {\tilde{f}}_{j} | = O (1 / \sqrt{N})$ , then

| E [g_{i} (y (k), n_{qy}) e (k)] | = O (\sqrt{logN / N}), n_{qy} > n_{c} + n_{\hat{f}} and i = 1, 2, \dots n_{a}^{'}

Proof

In order to prove the assertion of the lemma, we first prove that $(E [(v (k) - \hat{v} (k)) 2]) 1 / 2$ is of order $O (\sqrt{logN / N})$ . We have the following decomposition

v (k) - \hat{v} (k) = (v (k) - \tilde{v} (k)) + (\tilde{v} (k) - \hat{v} (k))

(21)

Therefore, by Minkowski inequality

(E [(v (k) - \hat{v} (k)) 2]) \frac{1}{2} \leq (E [(v (k) - \tilde{v} (k)) 2]) \frac{1}{2} + (E [(\tilde{v} (k) - \hat{v} (k)) 2]) \frac{1}{2} {= I}_{k} + {II}_{k} (say)

(22)

Since $ɛ (k)$ be stationary ergodic with $E [ɛ (k)] = 0$ and $E [ɛ 2 (k)] = σ_{ɛ}^{2} < \infty$ . Therefore

I_{k} = (E [(\sum_{j = 1}^{\infty} f_{j} ɛ (k - j) - \sum_{j = 1}^{\infty} {\tilde{f}}_{j} ɛ (k - j)) 2]) \frac{1}{2} = (E [(\sum_{j = 1}^{\infty} (f_{j} - {\tilde{f}}_{j}) ɛ (k - j)) 2]) \frac{1}{2}

Again using Minkowski inequality

I_{k} \leq \sum_{j = 1}^{\infty} (E [((f_{j} - {\tilde{f}}_{j}) ɛ (k - j)) 2]) \frac{1}{2} = \sum_{j = 1}^{\infty} | f_{j} - {\tilde{f}}_{j} | (E [ɛ 2 (k - j)]) \frac{1}{2} = \sum_{j = 1}^{\infty} | f_{j} - {\tilde{f}}_{j} | σ_{ɛ}

(23)

Since $\sum {no}_{j = 0}^{\infty} j β | f_{j} | < \infty$ and $\sum {no}_{j = 0}^{\infty} j β | {\tilde{f}}_{j} | < \infty$ , therefore $\sum {no}_{j = 0}^{\infty} j β | {\tilde{f}}_{j} - f_{j} | < \infty$ . Hence given a positive number δ (however small), ∃ $τ_{1} \in ℤ$ ⁺ such that $\sum {no}_{j = τ_{1} + 1}^{\infty} | f_{j} - {\tilde{f}}_{j} | < \frac{δ}{σ_{ɛ}}$ . But then

I_{k} \leq (\sum_{j = 1}^{τ_{1}} | f_{j} - {\tilde{f}}_{j} | σ_{ɛ} + \sum_{j = τ_{1} + 1}^{\infty} | f_{j} - {\tilde{f}}_{j} | σ_{ɛ}) \leq {sup}_{j \in ℤ +} | {\tilde{f}}_{j} - f_{j} | (τ_{1} σ_{ɛ}) + δ = O (\sqrt{log N / N})

(24)

Again, by construction $| {\tilde{f}}_{j} - {\hat{f}}_{j} | = 0, \forall j \leq n_{\hat{f}}$ , therefore

{II}_{k} = (E [(\sum_{j = 1}^{\infty} {\tilde{f}}_{j} ɛ (k - j) - \sum_{j = 1}^{n_{\hat{f}}} {\hat{f}}_{j} ɛ (k - j)) 2]) \frac{1}{2} = (E [(\sum_{j = n_{\hat{f}} + 1}^{\infty} {\tilde{f}}_{j} ɛ (k - j)) 2]) \frac{1}{2} \leq \sum_{j = n_{\hat{f}} + 1}^{\infty} (E [({\tilde{f}}_{j} ɛ (k - j)) 2]) \frac{1}{2} = \sum_{j = n_{\hat{f}} + 1}^{\infty} | {\tilde{f}}_{j} | σ_{ɛ} = O (1 / \sqrt{N})

(25)

Here we have again used the fact that ∃ $τ_{2} \in ℤ$ ⁺ such that $\sum {no}_{j = τ_{2} + 1}^{\infty} | {\tilde{f}}_{j} | < \frac{δ}{σ_{ɛ}}$ . Combining equations (22), (24), and (25) and using $(1 / \sqrt{N}) < (\sqrt{log N / N})$ we get

(E [(v (k) - \hat{v} (k)) 2]) 1 / 2 = O (\sqrt{log N / N}) + O (1 / \sqrt{N}) = O (\sqrt{log N / N})

(26)

Let $r_{ɛ v} (- t) = E [ɛ (k) v (k + t)]$ and $r_{ɛ \hat{v}} (- t) = E [ɛ (k) \hat{v} (k + t)]$ . By Cauchy–Schwartz inequality

| r_{ɛ v} (- t) - r_{ɛ \hat{v}} (- t) | \leq (E [ɛ 2 (k)]) \frac{1}{2} (E [(v (k + t) - \hat{v} (k + t)) 2]) \frac{1}{2}

(27)

Using equation (26), $| r_{ɛ v} (- t) - r_{ɛ \hat{v}} (- t) | = O (\sqrt{log N / N})$ for all positive t. Note that in equation (20), $ɛ (k)$ have finite dependence on $\hat{v} (k)$ and therefore $r_{ɛ \hat{v}} (- t) = 0$ ∀ $t > n_{\hat{f}}$ . Hence, for $t > n_{\hat{f}}$ , $| r_{ɛ v} (- t) | = O (\sqrt{log N / N})$ , Similarly $| r_{vv} (- t) | = O (\sqrt{log N / N})$ . Using equation (17)

| r_{ee} (- t) | = | \sum_{j_{1} = 0}^{n_{c}} \sum_{j_{2} = 0}^{n_{c}} c_{j_{1}} c_{j_{2}} E [v (k - j_{1}) v (k + t - j_{2})] | : (c_{0} = 1) = | \sum_{j_{1} = 0}^{n_{c}} \sum_{j_{2} = 0}^{n_{c}} c_{j_{1}} c_{j_{2}} E [v (k - j_{1}) v (k - j_{1} + t - j_{2} + j_{1})] | = | \sum_{j_{1} = 0}^{n_{c}} \sum_{j_{2} = 0}^{n_{c}} c_{j_{1}} c_{j_{2}} r_{vv} (- (t - j_{2} + j_{1})) |

Now $| r_{ee} (- t) | = O (\sqrt{log N / N})$ for $t > n_{c} + n_{\hat{f}}$ and using equation (4), the proof follows. ▪

4.3. Hyperbolic tangent nonlinear transformation and adequacy of the fit

The hyperbolic tangent function is considered for the nonlinear transformation such that

E [tanh (y (k - i - n_{qy})) e (k)] < E [y (k - i - n_{qy}) e (k)]

for some

n_{qy} > 0

(where

i = 1, 2, \dots, n_{a}'

). Since the integral involved in

E [tanh (y (k - i - n_{qy})) e (k)]

is computationally difficult, its numerical value could be bounded by the inequality derived as

E [(μ \tanh (y (k - i - n_{qy})) - e (k)) 2] \geq 0 for μ > 0

(28)

Since the left-hand side of inequality is nonnegative for all $μ > 0$ , therefore when equated to zero, the quadratic equation either has both roots are real and equal, or both roots are imaginary. Thus, the discriminant is always non positive which leads to

| E [\tanh (y (k - i - n_{qy})) e (k)] | \leq (E [tanh 2 (y (k - i - n_{qy}))]) 1 / 2 (E [e 2 (k)]) 1 / 2

(29)

Now the expected value $E [tanh 2 (y (k - i - n_{qy}))]$ and $E [e 2 (k)]$ can be replaced with $\frac{1}{N} \sum {no}_{k = 1}^{N} tanh 2 (y (k - i - n_{qy})$ and $\frac{1}{N} \sum {no}_{k = 1}^{N} e 2 (k)$ (Terrence and Andrew, 2010). Since, $| tanh x | \leq 1 \forall x \in ℝ$ , it follows that $0 \leq \frac{1}{N} \sum {no}_{k = 1}^{N} tanh 2 (y (k - i - n_{qy}) \leq 1$ for all N and therefore also holds when $N \to \infty$ . Hence

| E [\tanh (y (k - i - n_{qy})) e (k)] | \leq (E [e 2 (k)]) 1 / 2

(30)

Thus, by using $tanh$ transformation we get a sharper bound.

An examination of the adequacy of the fit is needed to give the way in which the extended dimension $n_{a}'$ and delay of outputs $n_{qy}$ is chosen (see equation (14)). Among different pairs $n_{qy}, n_{a}'$ of positive integer, where, $n_{a}' \geq n_{a}$ , a pair can be searched such that $E [q (k) e (k)] = 0$ and the corresponding estimate $\hat{θ} IV$ provides the best fit. Then

r_{q ρ} = E [q (k) (y (k) - z T (k) \hat{θ} IV)] = E [q (k) (z T (k) θ \circ + e (k) - z T (k) \hat{θ} IV)] = E [q (k) z T (k)] (θ \circ - \hat{θ} IV) \to 0 as \hat{θ} IV \to θ \circ

(31)

For the best solution $\hat{θ} IV$ is chosen such that the sample cross product between instrumental vector and residual is close to zero, i.e. ${\hat{r}}_{q ρ} \to 0$ (see equation (9b)). This amounts to first test the hypothesis that $E [q (k) e (k)] = 0$ . For this fix $n_{a}' = n_{a}$ and increase the value of $n_{qy}$ and stop when minimum is achieved. After this, increase $n_{a}'$ to find the efficient estimates.

5. Analysis and consistency properties

The EIV method using nonlinear instruments does not involve the prefiltering operation or the augmented model, and thus does not require iterative procedures. Therefore, the consistency can be shown by conditions of equations (9) and (11). Using lemma 1 and lemma 2, the condition of equation (9) is already verified. In order to prove equation (11), partition the matrix ${\hat{R}}_{qz}$ as

{\hat{R}}_{qz} = \frac{1}{N} [Q_{y} Q_{f (u)}] [Z_{y} Z_{f (u)}] = \frac{1}{N} [Q_{y} Z_{y} Q_{y} Z_{f (u)} Q_{f (u)} Z_{y} Q_{f (u)} Z_{f (u)}] Q_{y} = [q_{y} (1), q_{y} (2), \dots q_{y} (N)], Q_{f (u)} = [q_{f (u)} (1), q_{f (u)} (2), \dots q_{f (u)} (N)] Z_{y} = [z_{y} (1), z_{y} (2) \dots, z_{y} (N)] T, Z_{f (u)} = [z_{f (u)} (1), z_{f (u)} (2), \dots z_{f (u)} (N)] T

(32)

where

q_{y} (k) = [g_{1} (y (k), n_{qy}) g_{2} (y (k), n_{qy}) \dots g_{n_{a}^{'}} (y (k), n_{qy})] T

and

z_{y} (k) = [- y (k - 1) - y (k - 2) \dots - y (k - n_{a})] T

. Also

q_{f (u)} (k) = z_{f (u)} (k) = [u (k - 1), \dots u n_{p} (k - 1) \dots \dots u (K - n_{b}), \dots u n_{p} (K - n_{b})] T

The $lim {NO}_{N \to \infty} {\hat{R}}_{qz}$ may or may not equal to $R_{qz}$ because the expectations are not preserved under nonlinear transformations. The explicit computation of $R_{qz}$ is difficult. In such cases, the matrix ${\hat{R}}_{qz}$ can be used to compute the approximate solution (Laurain et al., 2008).

Theorem 1

Let M be an idempotent matrix such that $rank (Q_{y} M Z_{y}) = n_{a}$ and $rows (Q_{y}) \geq coluoumns (Z_{y})$ . If input is strongly nonlinearly persistently excited, then the estimates provided by the IV methods will satisfy the consistency condition of equation (11).

Proof

To prove that $rank ({\hat{R}}_{qz}^{T} {\hat{R}}_{qz}) = n_{a} + n_{b} n_{p},$ it is sufficient to prove that $rank ({\hat{R}}_{qz}) = n_{a} + n_{b} n_{p}$ (Albert, 1972). Now pre-multiplying and post-multiplying the matrix ${\hat{R}}_{qz}$ by nonsingular matrices, we get

\frac{1}{N} [I - Q_{y} Z_{f (u)} (Q_{f (u)} Z_{f (u)}) - 1 0 I] [Q_{y} Z_{y} Q_{y} Z_{f (u)} Q_{f (u)} Z_{y} Q_{f (u)} Z_{f (u)}] \times [I 0 - (Q_{f (u)} Z_{f (u)}) - 1 Q_{y} Z_{y} I] = \frac{1}{N} [Q_{y} Z_{y}^{'} - Q_{y} Z_{f (u)} (Q_{f (u)} Z_{f (u)}) - 1 Q_{f (u)} Z_{y} 00 Q_{f (u)} Z_{f (u)}]

This implies

\frac{1}{N} [Q_{y} Z_{y} Q_{y} Z_{f (u)} Q_{f (u)} Z_{y}^{'} Q_{f (u)} Z_{f (u)}] = \frac{1}{N} [I Q_{y} Z_{f (u)} (Q_{f (u)} Z_{f (u)}) - 1 0 I] [A' 00 Q_{f (u)} Z_{f (u)}] \times [I 0 (Q_{f (u)} Z_{f (u)}) - 1 Q_{y} Z_{y} I]

where

A' = Q_{y} Z_{y} - Q_{y} Z_{f (u)} (Q_{f (u)} Z_{f (u)}) - 1 Q_{f (u)} Z_{y}

. Since pre-multiplication and post-multiplication by a nonsingular matrix does not alter the rank of the matrix, therefore

rank [Q_{y} Z_{y} Q_{y} Z_{f (u)} Q_{f (u)} Z_{y} Q_{f (u)} Z_{f (u)}] = rank (A') + rank (Q_{f (u)} Z_{f (u)})

The construction of instrumental vector given by (14) implies that $rank (Q_{f (u)} Z_{f (u)}) =$ $rank (Z_{f (u)}^{T} Z_{f (u)})$ . Now if the input is strongly nonlinearly persistently excited, then $Z_{f (u)}^{T} Z_{f (u)} ≻ 0$ . So $rank (Q_{f (u)} Z_{f (u)}) = n_{b} n_{p}$ and

rank (A') = rank (Q_{y} (I - Z_{f (u)} (Z_{f (u)}^{T} Z_{f (u)}) - 1 Z_{f (u)}^{T}) Z_{y})

Since the matrix $Z_{f (u)}^{T} Z_{f (u)} ≻ 0$ therefore, $rank (A') = rank (Q_{y} M Z_{y})$ , where $M = I - Z_{f (u)} (Z_{f (u)}^{T} Z_{f (u)}) - 1 Z_{f (u)}^{T}$ is an idempotent matrix. Hence if $rank (Q_{y} M Z_{y}) = n_{a}$ , then $rank ({\hat{R}}_{qz}) = n_{a} + n_{b} n_{p}$ thereby proving the condition of equation (11). ▪

Next it is assumed that $lim {no}_{N \to \infty} {\hat{R}}_{qz}^{T} {\hat{R}}_{qz}$ exists and is nonsingular, and $Γ = lim {no}_{N \to \infty} ({\hat{R}}_{qz}^{T} {\hat{R}}_{qz}) - 1 {\hat{R}}_{qz}^{T}$ , then the asymptotic efficiency of the estimates is given by:

Theorem 2

Suppose assumptions $A_{1}$ and $A_{2}$ hold, and $lim {no}_{N \to \infty} {\hat{R}}_{qz}^{T} {\hat{R}}_{qz}$ is nonsingular. If the instrumental vector $q T (k)$ defined by equation (14) is independent of the noise $e (k)$ , then the IV estimate is asymptotically normally distributed with zero mean and covariance matrix P, i.e.

\sqrt{N} (\hat{θ} IV - θ \circ) \to dist N ormal (0, P)

where

P = Γ (\sum_{τ = - n_{c}}^{n_{c}} E [q (k + τ) q T (k)] E [e (k) e (k + τ)]) Γ T, when D (q - 1) = 1 P = Γ (\sum_{τ = - (n_{c} + n_{\hat{f}})}^{(n_{c} + n_{\hat{f}})} E [q (k + τ) q T (k)] E [e (k) e (k + τ)]) Γ T, when D (q - 1) \neq 1

Proof

The parameter estimation error is given by

\hat{θ} IV - θ \circ = ({\hat{R}}_{qz}^{T} {\hat{R}}_{qz}) - 1 {\hat{R}}_{qz}^{T} [\frac{1}{N} \sum_{k = 1}^{N} q (k) y (k) - \frac{1}{N} \sum_{k = 1}^{N} q (k) z T (k) θ \circ] = ({\hat{R}}_{qz}^{T} {\hat{R}}_{qz}) - 1 {\hat{R}}_{qz}^{T} {\hat{r}}_{qe}

Now to prove the result, it is sufficient to show that for any positive integers k and s

{lim}_{N \to \infty} \frac{1}{N} E [\sum_{k = 1}^{N} q (k) e (k) \sum_{s = 1}^{N} q T (s) e (s)] = \sum_{τ = - n_{c}}^{n_{c}} E [q (k + τ) q T (k)] E [e (k) e (k + τ)]

By assumptions A₁, A₂, and lemma 1, $e (k)$ is independent of $q (k), q (s)$ and $e (s)$ for $k > s + n_{c}$ , therefore the above equality holds (Soderstrom and Stoica, 1981). Similarly, for model structure $D (q - 1) \neq 1$ , assumptions A₁, A₂ and lemma 2 shows that the equality holds for $k > s + n_{c} + n_{\hat{f}}$ with a slight modification wherein the right-hand side is given by

\sum_{τ = - (n_{c} + n_{\hat{f}})}^{(n_{c} + n_{\hat{f}})} E [q (k + τ) q T (k)] E [e (k) e (k + τ)] + O (log N / N)

Now the proof follows using the convergence results given by (Caines, 1976).

6. Simulation examples

In this section, two simulation examples are presented to illustrate the merits of nonlinear transformed instruments for the Hammerstein model. The input $u (k)$ is a stochastic sequence having zero mean unit variance with Normal distribution covering the sufficient range of the nonlinearity. Each example has been investigated with 3000 input–output measurements and two different noise levels corresponding to different signal to noise ratio (SNR). The stochastic properties are investigated using 300 Monte Carlo simulations with different realization of input $u (k)$ and noise $η (k)$ . The accuracy and convergence analysis are evaluated based on the performance criterion i.e. statistically estimation error (SEE) and average coefficient of variation (ACV). The SEE is defined as

SEE = \sqrt{\frac{m ({\hat{θ}}_{s}) - θ_{s} \circ 2}{θ_{s} \circ 2}}

(33)

where

m ({\hat{θ}}_{s})

represents the sample mean of estimator for 300 simulation runs. And, the ACV is computed as (

{\hat{θ}}_{s} (i)

represents the i th component)

ACV = \frac{1}{n_{s}} \sum_{i = 1}^{n_{s}} \frac{σ ({\hat{θ}}_{s} (i))}{| m ({\hat{θ}}_{s} (i)) |}

(34)

where

n_{s} = n_{a} + n_{b} + n_{p} - 1

and

σ ({\hat{θ}}_{s} (i))

represents the standard deviation of estimates.

6.1. Example 1

The discrete Hammerstein ARMAX model to be identified is given by

y (k) = \frac{0.80 q - 1 - 0.60 q - 2 + 0.95 q - 3}{1 - 1.6 q - 1 + 1.3 q - 2 - 0.28 q - 3} f (u (k)) + η (k)

where

f (u (k)) = 1 u (k) + 1.5 u 2 (k) + 0.4 u 3 (k)

and the noise

η (k)

is described by

η (k) = \frac{1 - 0.75 q - 1 + 0.90 q - 2 - 0.45 q - 3}{1 - 1.6 q - 1 + 1.3 q - 2 - 0.28 q - 3} ɛ (k)

In this example, the proposed Hammerstein nonlinear extended IV method (HNEIV) method is compared with the LS method, the nonlinear LS method (NLS), the ELS method, the DFLS technique based on the bilinear parameterization (DFLSBP), the iterative DFLS method based on the partial substitution (DFLSPS), and the IV method with linear delayed outputs briefly written as Hammerstein linear EIV (HLEIV). The DFLSBP and the DFLSPS algorithms are implemented using the ARMAX structure with MA noise model (i.e.

n_{c} = 3

). For iterative DFLSBP algorithm we closely follow the steps given in Wang and Ding (2016a) and Wang and Ding (2015). The DFLSPS method is implemented based on the partial substitution of the nonlinearity

f (u (k))

using the b₁ term (Wang and Ding, 2015). Further, the NLS method is implemented using “lsqnonlin” function of the Matlab optimization toolbox. The statistical properties for

SNR =

14.31 dB and 06.99 dB are presented in Table 1 and Table 2, respectively.

Table 1.

Comparative performance of parameter estimates for example 1 ( $N = 3000, SNR = 14.31 dB, 300$ independent runs $)$ .

Parameters		a ₁	a ₂	a ₃	b ₁	b ₂	b ₃	p ₂	p ₃	SEE	ACV
LS	mean	–1.5298	1.1925	–0.2142	0.7996	–0.5459	0.9111	1.5038	0.4008	0.2325	0.0315
LS	std	0.0086	0.0134	0.0083	0.0272	0.0276	0.0267	0.0503	0.0198	0.2325	0.0315
NLS	mean	–1.5316	1.1950	–0.2163	0.8031	–0.5361	0.9155	1.5164	0.4008	0.2327	0.0273
NLS	std	0.0085	0.0131	0.0081	0.0232	0.0177	0.0265	0.0447	0.0177	0.2327	0.0273
ELS	mean	–1.5847	1.2772	–0.2645	0.8009	–0.5898	0.9421	1.5025	0.4000	0.1079	0.0245
ELS	std	0.0054	0.0080	0.0053	0.0205	0.0240	0.0268	0.0433	0.0171	0.1079	0.0245
DFLSPS	mean	–1.6016	1.3032	–0.2811	0.7982	–0.5993	0.9482	1.5095	0.4031	0.0612	0.0391
DFLSPS	std	0.0039	0.0057	0.0038	0.0417	0.0312	0.0494	0.0817	0.0331	0.0612	0.0391
DFLSBP	mean	–1.6004	1.3008	–0.2803	0.7992	–0.6003	0.9487	1.5022	0.4010	0.0321	0.0153
DFLSBP	std	0.0034	0.0050	0.0033	0.0160	0.0137	0.0162	0.0269	0.0106	0.0321	0.0153
HLEIV	mean	–1.6029	1.3033	–0.2822	0.8000	–0.6043	0.9509	1.4997	0.4005	0.0475	0.0429
HLEIV	std	0.0310	0.0365	0.0241	0.0281	0.0371	0.0291	0.0506	0.0198	0.0475	0.0429
HNEIV	mean	–1.5996	1.2995	–0.2796	0.8001	–0.6015	0.9504	1.5000	0.4005	0.0248	0.0384
HNEIV	std	0.0233	0.0271	0.0189	0.0281	0.0340	0.0281	0.0505	0.0197	0.0248	0.0384
True value		–1.60	1.3	–0.28	0.80	–0.60	0.95	1.50	0.4

SNR: signal to noise ratio; SEE: statistically estimation error; ACV: average coefficient of variation; LS: least square; NLS: nonlinear least square; ELS: extended least square; DFLSPS: data filtering-based least square method based on the partial substitution; DFLSBP: data filtering-based least square technique based on the bilinear parameterization; HLEIV: Hammerstein linear extended instrumental variable; HNEIV: Hammerstein nonlinear extended IV method.

Table 2.

Comparative performance of parameters estimates for example 1 ( $N = 3000, SNR = 06.99 dB, 300$ independent runs $)$ .

Parameters		a ₁	a ₂	a ₃	b ₁	b ₂	b ₃	p ₂	p ₃	SEE	ACV
LS	mean	–1.3351	0.8983	–0.0294	0.8004	–0.3900	0.8108	1.5166	0.4061	0.4512	0.1614
LS	std	0.0219	0.0331	0.0201	0.0611	0.0649	0.0662	0.1362	0.0568	0.4512	0.1614
NLS	mean	–1.3388	0.9034	–0.0340	0.8093	–0.3579	0.8185	1.5854	0.4015	0.4547	0.1323
NLS	std	0.0219	0.0331	0.0200	0.0562	0.0338	0.0582	0.1138	0.0444	0.4547	0.1323
ELS	mean	–1.5268	1.1947	–0.2086	0.8015	–0.5439	0.9162	1.5229	0.4025	0.2350	0.0647
ELS	std	0.0146	0.0212	0.0132	0.0517	0.0577	0.0671	0.1106	0.0444	0.2350	0.0647
DFLSPS	mean	–1.5532	1.2305	–0.2287	0.7921	–0.5649	0.9362	1.5595	0.4049	0.2031	0.0734
DFLSPS	std	0.0145	0.0219	0.0139	0.0701	0.0508	0.0822	0.1423	0.0573	0.2031	0.0734
DFLSBP	mean	–1.5811	1.2888	–0.2659	0.7928	–0.6063	0.9461	1.5280	0.4026	0.1163	0.0341
DFLSBP	std	0.0069	0.0103	0.0069	0.0373	0.0292	0.0385	0.0588	0.0241	0.1163	0.0341
HLEIV	mean	–1.6072	1.3081	–0.2854	0.8008	–0.6072	0.9503	1.5138	0.4044	0.0828	0.1127
HLEIV	std	0.0844	0.1014	0.0682	0.0688	0.0948	0.0737	0.1282	0.0518	0.0828	0.1127
HNEIV	mean	–1.5975	1.2961	–0.2776	0.8006	–0.5993	0.9479	1.5055	0.4043	0.0554	0.0806
HNEIV	std	0.0326	0.0413	0.0284	0.0665	0.0725	0.0699	0.1282	0.0516	0.0554	0.0806
True value		–1.60	1.3	–0.28	0.80	–0.60	0.95	1.50	0.4

From Tables 1 and 2, it is observed that the estimates obtained from the LS and NLS methods are inconsistent and heavily biased. In both SNR cases, the ELS estimates are biased which is evident from the SEE values. From the results, it follows that the DFLSPS and DFLSBP algorithms provides consistent estimates with lower ACV in the low noise case (SNR = $1$ 4.31 dB) than the other methods. However, for SNR = 06.99 dB, both algorithms lose the unbiasedness property despite of lower ACV. The HLEIV method provides reasonable accuracy and manages to reduce the SEE which is still valuable in low SNR. On the contrary, HNEIV provides more accurate estimates with minimum SEE for both noise cases.

6.2. Example 2

In this example, the entire configuration is the same as for example 1 except that the noise is generated by a low-pass IIR digital filter. The objective of this example is to illustrate the robustness of nonlinear instruments when the structure obtained by combining stochastic noise part with Hammerstein system does not belong to the ARMAX model structure. The deterministic part i.e. Hammerstein system without noise has no structural errors. In particular, the noise $η (k)$ is given by

η (k) = H (z - 1) ɛ (k)

where

H (z - 1)

is the low pass IIR filter of 10th order with cut-off frequency equal to 0.12.

Monte Carlo simulations of 300 independent runs have been performed with two different SNR values, 14.22 dB and 07.07 dB. The noise model order for the ELS method is taken as 7th and 8th order polynomial, i.e. $n_{c} = 7, 8$ . The DFLSBP method is implemented using two different noise model structures namely the MA noise ( $n_{c} > 0, n_{d} = 0$ ) and the ARMA noise ( $n_{c} > 0, n_{d} > 0$ ). The choice of noise model order is selected among different orders based on the best performance criterion i.e. one that has minimum SEE. The MA noise model orders for both noise levels using the DFLSBP method are found as $n_{c} = 7$ . However, when the DFLSBP algorithm is implemented using the ARMA noise model, the model orders are obtained as $n_{c} = 2, n_{d} = 2$ for SNR = 14.22 dB and $n_{c} = 3, n_{d} = 2$ for SNR = 07.07 dB.

For HLEIV and HNEIV, two groups of estimates of parameters are obtained using two different combinations of delay in the outputs ( $n_{qy}$ ) and number of instruments ( $n_{a}'$ ), that has best statistical properties. For this purpose, the parameter estimates are calculated using combinations of $n_{qy}$ and $n_{a}'$ ranges between 0 and 30. The SEE values obtained from different combinations for SNR = 14.22 dB and 07.07 dB using the HLEIV and the HNEIV methods are illustrated in Figures 2 and 3, respectively. The statstical performance for $SNR =$ 14.22 dB and 07.07 dB using different approaches are listed in Tables 3 and 4, respectively. Further, the box distributions of the estimates using the HNEIV method are illustrated in Figure 4. Note that the SEE values denoted by ‘Error’ represented in Figures 2 and 3 are evaluated using the estimates ( $\hat{θ}$ ) of augmented parameter vector in equation (7b).

Table 3.

Comparative performance of parameters estimates for example 2 ( $N = 3000, SNR = 14.22 dB, 300$ independent runs $)$ .

Parameters		a ₁	a ₂	a ₃	b ₁	b ₂	b ₃	p ₂	p ₃	SEE	ACV
LS	mean	–1.6600	1.3678	–0.3294	0.7994	–0.6485	0.9630	1.4951	0.4004	0.1976	0.0149
LS	std	0.0077	0.0094	0.0071	0.0189	0.0179	0.0196	0.0081	0.0035	0.1976	0.0149
ELS $n_{c} = 7$	mean	–1.6238	1.3155	–0.3000	0.8001	–0.6194	0.9459	1.4995	0.4003	0.1167	0.0131
ELS $n_{c} = 7$	std	0.0053	0.0075	0.0060	0.0138	0.0171	0.0207	0.0054	0.0023	0.1167	0.0131
ELS $n_{c} = 8$	mean	–1.6234	1.3157	–0.3005	0.7998	–0.6191	0.9463	1.4994	0.4003	0.1165	0.0130
ELS $n_{c} = 8$	std	0.0056	0.0079	0.0060	0.0136	0.0165	0.0207	0.0054	0.0023	0.1165	0.0130
DFLSBP $n_{c} = 7$	mean	–1.6322	1.3310	–0.3049	0.7989	–0.6249	0.9637	1.5006	0.4003	0.1411	0.0158
DFLSBP $n_{c} = 7$	std	0.0115	0.0161	0.0106	0.0134	0.0201	0.0147	0.0054	0.0019	0.1411	0.0158
DFLSBP $(n_{c} = 2,$ $n_{d} = 2)$	mean	–1.6050	1.3045	–0.2843	0.8013	–0.6012	0.9544	1.5000	0.4000	0.0562	0.0019
DFLSBP $(n_{c} = 2,$ $n_{d} = 2)$	std	0.0014	0.0015	0.0014	0.0017	0.0022	0.0023	0.0002	0.0001	0.0562	0.0019
HLEIV $(n_{qy} = 8,$ $n_{a}' = 4)$	mean	–1.5998	1.3058	–0.2840	0.7985	–0.6008	0.9538	1.5008	0.4007	0.0530	0.0595
HLEIV $(n_{qy} = 8,$ $n_{a}' = 4)$	std	0.0649	0.0832	0.0578	0.0215	0.0552	0.0313	0.0098	0.0040	0.0530	0.0595
HLEIV $(n_{qy} = 7,$ $n_{a}' = 7)$	mean	–1.5854	1.2776	–0.2595	0.7983	–0.5893	0.9416	1.5029	0.4007	0.1114	0.0548
HLEIV $(n_{qy} = 7,$ $n_{a}' = 7)$	std	0.0527	0.0687	0.0489	0.0232	0.0465	0.0301	0.0140	0.0055	0.1114	0.0548
HNEIV $(n_{qy} = 9,$ $n_{a}' = 17)$	mean	–1.6028	1.2992	–0.2776	0.7984	–0.6034	0.9463	1.5012	0.4006	0.0475	0.0481
HNEIV $(n_{qy} = 9,$ $n_{a}' = 17)$	std	0.0494	0.0650	0.0448	0.0219	0.0419	0.0302	0.0086	0.0031	0.0475	0.0481
HNEIV $(n_{qy} = 9,$ $n_{a}' = 17)$	mean	–1.6029	1.3042	–0.2758	0.7985	–0.6037	0.9502	1.5021	0.4006	0.0523	0.0456
HNEIV $(n_{qy} = 9,$ $n_{a}' = 17)$	std	0.0456	0.0588	0.0421	0.0222	0.0400	0.0293	0.0085	0.0031	0.0523	0.0456
True value		–1.60	1.3	–0.28	0.80	–0.60	0.95	1.50	0.4

Table 4.

Comparative performance of parameters estimates for example 2 ( $N = 3000, SNR = 07.07 dB, 300$ independent runs $)$ .

Parameters		a ₁	a ₂	a ₃	b ₁	b ₂	b ₃	p ₂	p ₃	SEE	ACV
LS	mean	–1.7991	1.5230	–0.4410	0.7977	–0.7601	0.9928	1.4794	0.4018	0.3588	0.0289
LS	std	0.0180	0.0209	0.0147	0.0374	0.0385	0.0407	0.0200	0.0089	0.3588	0.0289
ELS $n_{c} = 7$	mean	–1.6856	1.3614	–0.3541	0.8008	–0.6682	0.9411	1.4933	0.4009	0.2230	0.0261
ELS $n_{c} = 7$	std	0.0117	0.0184	0.0130	0.0264	0.0351	0.0403	0.0128	0.0056	0.2230	0.0261
ELS $n_{c} = 8$	mean	–1.6848	1.3601	–0.3537	0.8016	–0.6674	0.9406	1.4930	0.4008	0.2218	0.0260
ELS $n_{c} = 8$	std	0.0117	0.0181	0.0126	0.0273	0.0355	0.0396	0.0129	0.0056	0.2218	0.0260
DFLSBP $n_{c} = 7$	mean	–1.7044	1.3826	–0.3519	0.7960	–0.6811	0.9987	1.5044	0.4010	0.2465	0.0260
DFLSBP $n_{c} = 7$	std	0.0181	0.0261	0.0179	0.0244	0.0399	0.0267	0.0088	0.0024	0.2465	0.0260
DFLSBP $(n_{c} = 3,$ $n_{d} = 2)$	mean	–1.6344	1.3288	–0.3076	0.8065	–0.6122	0.9757	1.5003	0.4001	0.1431	0.0116
DFLSBP $(n_{c} = 3,$ $n_{d} = 2)$	std	0.0115	0.0137	0.0105	0.0071	0.0112	0.0110	0.0021	0.0007	0.1431	0.0116
HLEIV $(n_{qy} = 8,$ $n_{a}' = 4)$	mean	–1.6115	1.3102	–0.2627	0.8019	–0.6093	0.9516	1.5043	0.4020	0.0932	0.1106
HLEIV $(n_{qy} = 8,$ $n_{a}' = 4)$	std	0.1034	0.1351	0.0999	0.0498	0.0944	0.0621	0.0346	0.0129	0.0932	0.1106
HLEIV $(n_{qy} = 7,$ $n_{a}' = 7)$	mean	–1.5801	1.2978	–0.2969	0.8070	–0.5875	0.9556	1.5747	0.4443	0.1770	0.4253
HLEIV $(n_{qy} = 7,$ $n_{a}' = 7)$	std	0.2366	0.2958	0.2046	0.1129	0.2290	0.1313	0.8435	0.5034	0.1770	0.4253
HNEIV $(n_{qy} = 8,$ $n_{a}' = 14)$	mean	–1.6037	1.2899	–0.2766	0.8015	–0.6012	0.9436	1.5055	0.4032	0.0704	0.0981
HNEIV $(n_{qy} = 8,$ $n_{a}' = 14)$	std	0.0900	0.1210	0.0840	0.0579	0.0913	0.0673	0.0225	0.0083	0.0704	0.0981
HNEIV $(n_{qy} = 9,$ $n_{a}' = 17)$	mean	–1.6159	1.2937	–0.2587	0.8018	–0.6102	0.9387	1.5062	0.4029	0.1044	0.0843
HNEIV $(n_{qy} = 9,$ $n_{a}' = 17)$	std	0.0703	0.0935	0.0671	0.0553	0.0783	0.0624	0.0230	0.0082	0.1044	0.0843
True value		–1.60	1.3	–0.28	0.80	–0.60	0.95	1.50	0.4

It is observed from Tables 3 and 4, for the case of SNR $=$ 14.22 dB, both methods (the ELS method and the DFLSBP method using MA noise model) provide biased estimates which are more prominent when the SNR further reduces to 07.07 dB. This is expected since in the low SNR case, the system combined with IIR noise does not belong to the ARMAX model. Therefore, the biases in the estimates are due to the misfit of the noise model. The estimates obtained using the DFLSBP method, when it is implemented using the ARMA noise model are consistent and provides lower SEE and ACV for the SNR = 14.22 dB. However, the estimates are slightly biased when the SNR reduces to 07.07 dB. The HNEIV method provides consistent estimates and reduces the SEE, which is significantly lower than the values obtained from other methods.

7. Conclusion

An extended instrumental variable method using the instruments obtained from the nonlinear transformation of observed variables for the estimation of Hammerstein ARMAX models has been proposed and analyzed. It has been shown that the proposed method can provide consistent estimates, even in case of Hammerstein models with general noise model structures. The consistency of the estimates for the case of general noise model has been shown by computing the bounds on the expected relationship between instrumental variables and noise. A significant advantage of the method is that, it eliminates the noise induced bias without explicitly incorporating the noise model and filtering operation, and thus has simpler structure.

In addition, the detailed simulation study of the case when the noise is generated through IIR filter has been performed. The simulation results show that the nonlinear transformation is more prominent as compared with the linear transformation in terms of relatively lower mean around true values and low variance. It also has been demonstrated that the method can perform in low SNR environment without losing too much estimation accuracy and statistical properties of the estimates. Furthermore, the comparison with the LS variants indicates that the method has offered better performance with relatively lower SEE.

Footnotes

Declaration of Conflicting Interests

The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding

The author(s) received no financial support for the research, authorship, and/or publication of this article.

References

Albert

(1972) Regression and the Moore–Penrose Pseudoinverse, New York: Academic Press.

Bai

Chan

(2008) Identification of an additive nonlinear system and its applications in generalized Hammerstein models. Automatica 44(2): 430–436.

Balestrino

Landi

Ould-Zmirli

et al. (2001) Automatic nonlinear auto-tuning method for Hammerstein modeling of electrical drives. IEEE Transactions on Industrial Electronics 48(3): 645–655.

Berk

(1974) Consistent autoregressive spectral estimates. The Annals of Statistics 2(3): 489–502.

Box

GEP

Pierce

(1970) Distribution of residual autocorrelations in Autoregressive-integrated moving average time series models. Journal of the American Statistical Association 65(332): 1509–1526.

Bruls

Chou

Haverkamp

et al. (1999) Linear and nonlinear system identification using separable least squares. European Journal of Control 5(11): 116–128.

Buhlmann

(1995) Moving-average representation of autoregressive approximations. Stochastic Processes and their Applications 60: 331–342.

Caines

(1976) On the asymptotic normality of instrumental variable and least squares estimators. IEEE Transactions on Automatic Control 21(4): 598–600.

Ding

Chen

(2005) Identification of Hammerstein nonlinear ARMAX systems. Automatica 41(9): 1479–1489.

10.

Ding

Shi

Chen

(2007) Auxiliary model-based least-squares identification methods for Hammerstein output-error systems. Systems & Control Letters 56(5): 373–380.

11.

Grimmett

Stirzaker

(2001) Probability and Random Processes, Oxford: Oxford University Press.

12.

Ding

(2013) Multistage least squares based iterative estimation for feedback nonlinear systems with moving average noises using the hierarchical identification principal. Nonlinear Dynamics 73(1): 583–592.

13.

Jurado

(2006) A method for the identification of solid oxide fuel cells using a Hammerstein model. Journal of Power Sources 154(1): 145–152.

14.

Kuo

Hosein

Farmanborda

(2012) Nonlinear output feedback control of a flexible link using adaptive neural network: controller design. Journal of Vibration and Control 19(11): 1690–1708.

15.

Laurain V, Gilson M, Garnier H, et al. (2008) Refined instrumental variable methods for identification of Hammerstein continuous time Box–Jenkins model. In: Proceedings of the 47th IEEE conference on decision and control, Cancun, Mexico, 9–11 December 2008, pp. 1386–1391. IEEE. DOI: 10.1109/CDC.2008.4738853.

16.

Ding

Hua

(2014) Maximum likelihood Newton recursive and the Newton iterative estimation algorithms for Hammerstein CARAR systems. Nonlinear Dynamics 75(1): 235–245.

17.

Ljung

(1999) System Identification: Theory for the User, Upper Saddle River, NJ: Prentice-Hall.

18.

Narendra

Gallman

(1966) An iterative method for the identification of nonlinear systems using Hammerstein model. IEEE Transactions on Automatic Control 11(3): 546–550.

19.

Schoukens

Widanage

Godfrey

et al. (2007) Initial estimates for the dynamics of a Hammerstein system. Automatica 43(7): 1296–1301.

20.

Stone

(1948) The generalized Weierstrass approximation theorem. Mathematics Magazine 21(5): 237–254.

21.

Soderstrom

Stoica

(1981) Comparison of some instrumental variable methods-consistency and accuracy aspects. Automatica 17(1): 101–115.

22.

Soderstrom

Stoica

(1983) Instrumental Variable Methods for System Identification, Berlin: Springer.

23.

Soderstrom

Stoica

(1989) System identification, Eaglewood Cliffs, NJ: Prentice-Hall.

24.

Srinivasan

Rengaswamy

Narasimhan

et al. (2005) Control loop performance assessment. 2. Hammerstein model approach for stiction diagnosis. Industrial & Engineering Chemistry Research 44(17): 6719–6728.

25.

Stoica

Soderstrom

(1983) Optimal instrumental variable estimation and approximate implementations. IEEE Transactions on Automatic Control 28(7): 757–771.

26.

Terrence

Andrew

(2010) Uniform convergence of Vapnik–Chervonenkis classes under ergodic sampling. The Annals of Probability 38(4): 1345–1367.

27.

Vörös

(1999) Iterative algorithm for parameter identification of Hammerstein systems with two-segment nonlinearities. IEEE Transactions on Automatic Control 44(11): 2145–2149.

28.

Wang

(2016) Hierarchical parameter estimation for a class of MIMO Hammerstein systems based on the reframed models. Applied Mathematics Letters 57: 13–19.

29.

Wang

Ding

(2016a) Parameter estimation algorithms for multivariable Hammerstein CARMA systems. Information Sciences 355-356: 237–248.

30.

Wang

Ding

Chu

(2013) Data filtering based recursive least square algorithm for Hammerstein systems using the key term separation principal. Information Sciences 222: 203–212.

31.

Wang

Liu

Ding

(2015) Highly efficient identification methods for dual-rate Hammerstein systems. IEEE Transactions on Control System Technology 23(5): 1952–1960.

32.

Wang

Ding

(2015) Iterative estimation for a non-linear IIR filter with moving average noise by means of the data filtering technique. IMA Journal of Mathematics Control and Information. Epub ahead of print 24 December 2015. DOI: 10.1093/imamci/dnv067.

33.

Wang

Ding

(2016b) Novel data filtering based parameter identification for multiple-input multiple-output systems using the auxiliary model. Automatica 71: 308–313.

34.

Wang

Ding

(2016c) Recursive least squares algorithm and gradient algorithm for Hammerstein–Wiener systems using the data filtering. Nonlinear Dynamics 84(2): 1045–1053.

35.

Wang

(2014) Data-filtering-based iterative identification methods for nonlinear FIR-MA systems. Journal of Vibration and Control 20(14): 2193–2201.

36.

Westwick

Kearney

(2001) Separable least squares identification of nonlinear Hammerstein models: application to stretch reflex dynamics. Annals of Biomedical Engineering 29(8): 707–718.

37.

Xie

Yang

(2010) Gradient-based iterative identification for nonuniform sampling output error systems. Journal of Vibration and Control 17(3): 471–478.

38.

Young

(1970) An instrumental variable method for real time identification of a noisy process. Automatica 6(2): 271–287.