## 3 System Identification
### WANARX Model  
The approximation capability associated with wavelet functions is well known, so a wavelet network is selected as a nonlinearity estimator for the NARX model. To compact the notation, the regressors vector \(\mathbf{x}_{k}\in\mathcal{R}^{m}\) is introduced:  
\[\mathbf{x}_{k}=(y_{k-1},\ldots,y_{k-n_{a}},u_{k-n_{b}},\ldots,u_{k-n_{k}-n_{b} +1})^{T} \tag{13}\]  
where \(n_{k}\) is the input delay. By rewriting Equation (11) according to the WANARX structure, one obtains:  
\[y_{k}=f(\mathbf{x}_{k})+(\mathbf{x}_{k}-\bar{\mathbf{x}})^{T}\mathbf{P} \boldsymbol{\theta}+l+e_{k} \tag{14}\]  
where \(\bar{\mathbf{x}}\) is the regressor mean, \(\boldsymbol{\theta}\) the \(p\)-by-1 vector of linear coefficients, \(\mathbf{P}\) an \(m\)-by-\(p\) projection matrix, \(l\) the scalar output offset, and  
\[f(\mathbf{x}_{k})=W(\mathbf{x}_{k})+S(\mathbf{x}_{k})\]  
the nonlinear function constituted by dilated and translated wavelets and dilated and translated scaling functions (to improve the regularity of the estimator [23]). The nonlinear block \(f(\mathbf{x}_{k})\) and the linear block \((\mathbf{x}_{k}-\bar{\mathbf{x}})^{T}\mathbf{P}\boldsymbol{\theta}+l\) are combined together in a series-parallel structure to improve the estimation results.  
Wavelets are functions located both in time and frequency. A finite energy signal can be decomposed into different frequency components by the superposition of functions obtained through scaling and translating an initial function \(\psi\) known as the mother wavelet [23]. A radial function \(\psi_{w}\) depending on the squared magnitude of vector \(x\) is chosen as the mother wavelet function for the relationship described by \(W(\cdot)\):  
\[W(\mathbf{x})=\sum_{i=1}^{d_{w}}w_{i}\psi_{w}(b_{i}(\mathbf{x}-\bar{\mathbf{x }})^{T}\mathbf{Q}-\mathbf{C}_{i}) \tag{15}\]  
where  
\[\psi_{w}(X)=(dim(X)-XX^{T})e^{-XX^{T}/2} \tag{16}\]  
and \(\mathbf{Q}\) is a \(m\)-by-\(q\) projection matrix, \(m\geq q\), \(w_{1},\ldots,w_{d_{w}}\) the wavelet coefficients, \(b_{1},\ldots,b_{d_{w}}\) the wavelet dilations parameters, \(\mathbf{C}_{1},\ldots,\mathbf{C}_{d_{w}}\) the 1-by-\(q\) row vectors of wavelet translations, \(m\) the dimension of the regressors vector, \(\bar{\mathbf{x}}\) the mean of \(\mathbf{x}\) and \(d_{w}\) the total number of  
Figure 1: GP-NARX scheme.  
wavelets. For the relationship described by \(S(\cdot)\), a radial function \(\psi_{s}\) is chosen as the scaling mother function:  
\[S(\mathbf{x})=\sum_{i=1}^{d_{s}}s_{i}\psi_{s}(d_{i}(\mathbf{x}-\mathbf{\bar{x}}) ^{T}\mathbf{Q}-\mathbf{E}_{i}) \tag{17}\]  
where  
\[\psi_{s}(X)=e^{-XX^{T}/2} \tag{18}\]  
in which \(s_{1},\ldots,s_{d_{s}}\) are the scaling coefficients, \(d_{1},\ldots,d_{d_{s}}\) the scaling dilations parameters, \(\mathbf{E}_{1},\ldots,\mathbf{E}_{d_{s}}\) the 1-by-\(q\) row vectors of scaling translations, and \(d_{s}\) the total number of scaling functions.  
Structures in Equations (15) and (17) are very similar to the ones adopted by a neural network: a the wavelet network can be seen as a one-hidden-layer neural network with \(\psi\) as the activation function of the hidden layer [24]. A schematic representation is reported in Figure 2.