Orthogonal Floating Search Algorithms: From The Perspective of Nonlinear System Identification

TL;DR

This paper introduces an orthogonal floating search framework for nonlinear system structure selection, enhancing existing algorithms by integrating orthogonal space concepts and ERR metrics, and demonstrating improved performance on benchmark systems.

Contribution

It adapts floating search algorithms with orthogonal space and ERR metrics, eliminating nesting effects and improving nonlinear system identification.

Findings

Framework effectively identifies correct system structures.

Existing feature selection methods can be tailored for nonlinear systems.

Demonstrated improved performance on benchmark nonlinear systems.

Abstract

The present study proposes a new Orthogonal Floating Search framework for structure selection of nonlinear systems by adapting the existing floating search algorithms for feature selection. The proposed framework integrates the concept of orthogonal space and consequent Error-Reduction-Ratio (ERR) metric with the existing floating search algorithms. On the basis of this framework, three well-known feature selection algorithms have been adapted which include the classical Sequential Forward Floating Search (SFFS), Improved sequential Forward Floating Search (IFFS) and Oscillating Search (OS). This framework retains the simplicity of classical Orthogonal Forward Regression with ERR (OFR-ERR) and eliminates the nesting effect associated with OFR-ERR. The performance of the proposed framework has been demonstrated considering several benchmark non-linear systems. The results show that most of the existing feature selection methods can easily be tailored to identify the correct system structure of nonlinear systems.

Equations

$$J({\mathcal{X}}^{\star })=\underset{\mathcal{X}\subset {\mathcal{X}}_{model}}{\max }J(\mathcal{X})$$

$$y(k)$$

$$n$$

$$y(k)$$

$$\text{where, }{x}_1(k)$$

$$w_1=x_1,w_i=x_i-\sum _{m=1}^{i-1}\frac{⟨w_m{x}_i⟩}{w_m^T{w}_m}$$

$$y(k)$$

$$ERR(x_i)=\frac{g_i^2{w}_i^T{w}_i}{y^{T}y},\text{where, }{g}_i=\frac{w_i^{T}y}{w_i^T{w}_i},i\in [1,n]$$

$$J({\mathcal{X}}_k)=\sum _{i=1}^{k}ERR(x_i),\forall x_i\in {\mathcal{X}}_k$$

$$x^{MS}=\{x_i|J({\mathcal{X}}_k\cup x_i)=\underset{i}{\max }J({\mathcal{X}}_k\cup x_i),x_i\in {\mathcal{X}}_{model}\setminus {\mathcal{X}}_k\},\forall i\in [1,n]$$

$$x^{LS}=\{x_i|J({\mathcal{X}}_k\setminus x_i)=\underset{i}{\max }J({\mathcal{X}}_k\setminus x_i),x_i\in {\mathcal{X}}_k\},\forall i\in [1,n]$$

$${\mathcal{X}}_o^{MS}=\{{\mathcal{X}}_o^i|J({\mathcal{X}}_k\cup {\mathcal{X}}_o^i)=\underset{i}{\max }J({\mathcal{X}}_k\cup {\mathcal{X}}_o^i),{\mathcal{X}}_o^i\in {\mathcal{X}}_{model}\setminus {\mathcal{X}}_k\}$$

$${\mathcal{X}}_o^{LS}=\{{\mathcal{X}}_o^i|J({\mathcal{X}}_k\setminus {\mathcal{X}}_o^i)=\underset{i}{\max }J({\mathcal{X}}_k\setminus {\mathcal{X}}_o^i),{\mathcal{X}}_o^i\in {\mathcal{X}}_k\}$$

$${\mathcal{S}}_1:y(k)=0.5+0.5y(k-1)+0.8u(k-2)-0.05y(k-2{)}^2+u(k-1{)}^2+e(k)$$

$$\Omega =\{{\mathcal{X}}_{{\xi }_{min}},{\mathcal{X}}_{{\xi }_{min}+1},\dots ,{\mathcal{X}}_{{\xi }_{max}}\}$$

$${\tau }_{i,\xi }=\{\begin{array}{l}1,\text{if}{x}_i\in {\mathcal{X}}_{\xi }\\ 0,\text{otherwise}\end{array},\forall {\mathcal{X}}_{\xi }\in \Omega$$

$$S7:y(k)=$$

$$+0.7582u(k-1)u(k-2)+0.1750u(k-1)u(k-3)+0.0864u(k-2{)}^2$$

$$+0.4916u(k-2)u(k-3)+0.0711u(k-3{)}^2-0.0375y(k-1)-0.0598y(k-2)$$

$$-0.0370y(k-3)-0.0468y(k-4)-0.0476y(k-1{)}^2-0.0781y(k-1)y(k-2)$$

$$-0.0189y(k-1)y(k-3)-0.0626y(k-1)y(k-4)-0.0221y(k-2{)}^2$$

$$-0.0617y(k-2)y(k-3)-0.0378y(k-2)y(k-4)-0.0041y(k-3{)}^2$$

$$-0.0543y(k-3)y(k-4)-0.0603y(k-4{)}^2+e(k)$$

$${\mathcal{S}}_8:w(k)$$

$$y(k)$$

$$\text{with}u(k)$$

$$\ddot{y}(t)+2\zeta {\omega }_n\dot{y}(t)+{\omega }_n^{2}y(t)+{\omega }_n^2\epsilon y(t{)}^3-u(t)=0$$

$$y(k)=$$

$$+1.9792\times 10^{-6}u(k-2)-0.34289y(k-1{)}^3+0.22813y(k-2)y(k-1{)}^2$$

$$-1.079\times 10^{-7}u(k-1)y(k-1{)}^2-0.11847y(k-1)y(k-2{)}^2$$

$$AIC({\mathcal{X}}_{\xi })$$

$$BIC({\mathcal{X}}_{\xi })$$

$$FPE({\mathcal{X}}_{\xi })$$

$$LILC({\mathcal{X}}_{\xi })$$

Full text

1 Introduction

Over the last several decades, researchers from various fields of science and engineering have developed several methods to construct mathematical models from measured input-output data ( popularly known as system identification)  [1, 2]. Since most of the practical systems are nonlinear, the first step in the process of fitting the models is the choice of the model amongst various available models (e.g., Volterra, Wiener, Polynomial/Rational) which can effectively capture the dynamics in input-output data. In this study, the focus is on the identification of nonlinear systems represented by polynomial nonlinear auto-regressive with exogenous inputs (NARX) models [1]. Note that the identification of system is a two stage process which involves 1) Determination of the significant terms and 2) Estimation of the corresponding parameters. While the parameters can be estimated using least-squares based algorithms, the determination of the significant terms is a challenging task and often referred to as ‘structure selection’.

A close examination of the structure selection problem in nonlinear system identification reveals that this problem has significant similarities with the feature selection problem encountered in the pattern recognition. In essence, both problems belong to a broader class of combinatorial optimization problem which is referred here as the attribute selection problem. To understand this further, consider a superset of ‘$n$’ number of ‘attributes’, denoted by: $\mathcal{X}_{model}=\{x_{1},x_{2},\dots x_{n}\}$, where any attribute (say ‘$x_{i}$’), could represent either a feature (in pattern recognition) or a term of the NARX model (in system identification). The goal of the attribute selection is to determine a subset of significant attributes, $\mathcal{X}^{\star}\subset\mathcal{X}_{model}$, by maximizing a suitable criterion function,‘$J(\cdotp)$’, as follows:

[TABLE]

An exhaustive search to solve such combinatorial problem is often intractable even for a moderate number of attributes, ‘$n$’, as it requires the evaluation of $2^{n}$ subsets.

Over the years several approaches have been proposed to address the structure selection problem of NARX models, e.g., see [3, 4, 5, 6, 7, 8]. Among these, the Orthogonal Forward Regression (OFR) [9, 10, 4] is, arguably, computationally the most effective. Perhaps, for this reason, it has been extensively applied in many applications [11, 1]. The central idea of this approach is to decouple the terms in the orthogonal space so that the contribution of each term to the output variance can individually be determined. From the search perspective, the OFR is a sequential greedy approach where the term with the highest performance metric, referred to as Error-Reduction-Ratio (ERR), is included in each step and it has proven to be very effective in various applications [11, 3]. However, it has been shown that under certain conditions the OFR may yield sub-optimal term subsets [12, 13, 14, 15, 16]. The reason for such suboptimal performance is often attributed to the performance metric, ERR. For example, Mao and Billings [12] have shown that for the given term, the value of ERR is dependent on the orthogonalization path, i.e., the order in which the terms are orthogonalized. This may lead to the selection of spurious terms, when there is an information overlap among the non-orthogonal terms, especially in the earlier stages of the search [15].

To alleviate these issues, several approaches have been proposed to extend/modify the classical OFR algorithm. These can be categorized mainly into two approaches: 1) Two-step approach [12, 17, 15] and 2) Improved performance metric [13, 14, 16].

In the two-step approach, OFR is applied first to identify the initial term subset which is further refined by a secondary search procedure. For example, in [12], the initial term subset is refined by applying Genetic Algorithm (GA) to identify the optimal orthogonalization path. However, the explosive growth of the search space limits the application of this approach, e.g., for ‘$n$’ number of terms, the total number of possible orthogonalization paths becomes ‘$n!$’. In comparison, if GA (or any suitable algorithm) is applied directly to select term subset the search space reduces from $n!$ to $2^{n}$ [18]. The ‘iterative OFR’ (iOFR) algorithm [15] uses each term in the initial term subset as a pivot to identify new, and possibly better, term subsets in the secondary search. In [17], a backward elimination approach is proposed as the secondary search to remove the spurious terms contained in the initial term subset. Most of the two-step approach, e.g. [17], are usually effective provided all the system/significant terms are identified in the initial term subset.

The focus of the second approach is to replace the performance metric, ‘ERR’ [13, 14, 16]. For this purpose, several performance metrics have been proposed based on Simulation Error [13], Mutual Information [14] and Ultra-orthogonal least squares [16]. While these metrics improve the search performance of OFR in many cases, they often increase the computational complexity of OFR.

It is worth to emphasize that the aforementioned approaches to improve OFR are based on the prevailing belief that the sub-optimal performance of OFR is due only to its performance metric, ERR. However, in this study, it is argued that this may be in part due to the uni-directional sequential search nature of OFR where terms are progressively added into the term subset until pre-specified ERR threshold is reached. This leads to the ‘nesting effect’, i.e., once a term is included in the subset, it cannot be removed from the term subset. This imposes a hierarchy of terms, whereas in practice, all significant terms are equally important. Note that the nesting effect is not new and it has been well-known among the feature selection community since the early seventies [19, 20, 21]. Over the years, state-of-the-art ‘floating’ feature selection algorithms have specifically been developed to address the nesting effect [22, 23, 24, 25]. Since there already exist a class of floating feature selection algorithms and given that both structure selection and feature selection share a common search objective, the question arises: Can the feature selection algorithms be adapted to address the structure selection problem?

The present study aims to address this apropos question. In particular, the proposed Orthogonal Floating Search Framework integrates the key properties of the classical Sequential Forward Floating Search (SFFS) and its variants such as Improved Forward Floating Search (IFFS) and Oscillating Search (OS) with the Error-Reduction-Ratio (ERR) to determine the structure of nonlinear systems. In each step of these algorithms, after adding a ‘significant’ term (discussed in Section 3.1), the selected terms are re-examined to identify and remove the ‘least significant’ term (discussed in Section 3.1). Thus, the spurious terms can eventually be replaced by the system/significant terms during the search process. The search performance is therefore significantly improved without the need for a second search procedure or a complicated performance measure. The efficacy of the proposed approach is demonstrated on several benchmark non-linear systems. The worst case scenario (non-persistent excitation) and the identification of a discrete-time model of a continuous-time system are also considered. The results of this investigation convincingly demonstrate that the feature selection algorithms can indeed be tailored for the purpose of structure selection.

The rest of the article is organized as follows: The polynomial NARX model and the OFR approach are discussed briefly in Section 2. Next, the proposed Orthogonal Floating Search Framework is discussed in detail in Section 3. The framework of this investigation is described in Section 4. Finally, the results of this investigation are discussed at length in Section 5.

2 Preliminaries

In the following, the polynomial NARX model and the orthogonal regression approach are briefly discussed.

2.1 The Polynomial NARX Model

The NARX model represents a non-linear system as a function of recursive lagged input and output terms as follows:

[TABLE]

where $y(k)$ and $u(k)$ respectively represent the output and input at time intervals $k$, $n_{y}$ and $n_{u}$ are corresponding lags and $F^{n_{l}}\{\cdotp\}$ is some nonlinear function of degree $n_{l}$.

The total number of possible terms or model size ($n$) of the NARX model is given by,

[TABLE]

This model is essentially linear-in-parameters and can be expressed as:

[TABLE]

$p_{y},q_{u}\geq 0$; $1\leq p_{y}+q_{u}\leq n_{l}$; $1\leq n_{y_{j}}\leq n_{y}$;$1\leq n_{u_{k}}\leq n_{u}$; $n_{l}$ is the degree of polynomial expansion; ‘$\mathcal{N}$’ denotes the total number of data points.

2.2 Orthogonal Regression

In this study, the orthogonalization is used to determine the significance of a particular term, which will be discussed in Section 3.1. It is, therefore, pertinent to briefly discuss the concept of orthogonal regression introduced by Billings et al. [9, 10, 4]. In the orthogonal regression, as the name suggests, the original terms are replaced by the corresponding orthogonal terms. Given that the orthogonal terms are decoupled, the significance of each term can be determined independently in the orthogonal space. This could be achieved by a simple Gram-Schmidt orthogonalization:

[TABLE]

Consequently, the NARX model (3) can be represented in the orthogonal space as follows:

[TABLE]

Given that, in the orthogonal space, $w_{i}^{T}w_{j}=0$ holds for $i\neq j$, the significance of the $i^{th}$ term, $x_{i}$, can be easily be determined using the Error-Reduction-Ratio (ERR) metric as follows:

[TABLE]

3 Proposed Orthogonal Floating Search Framework

This study aims to investigate whether the existing feature selection algorithms can be tailored to address the non-linear system identification problem. For this purpose, a new search framework is being proposed which integrates the term orthogonalization into ‘floating’ feature selection algorithms and therefore referred to as ‘Orthogonal Floating Search Framework’.

The term orthogonalization decouples the structure selection and the parameter estimation, which are crucial steps of the nonlinear system identification. It is, therefore, a vital first step to adapt the feature selection algorithms for system identification. Further, the following ‘floating’ feature selection algorithms have been adapted: Sequential Forward Floating Search (SFFS) [22], Improved Forward Floating Search (IFFS) [25] and Oscillating Search (OS) [24]. The motivation for the selection of these methods is two-fold: 1) The term orthogonalization can be integrated with relative ease into the floating search methods. 2) The step-wise subset building approach of the sequential floating methods (e.g., SFFS and IFFS) is similar to that of the Orthogonal Forward Regression (OFR-ERR) [4]. Further, in OS, the depth of the search process can easily be controlled as per the prevailing requirements, as will be discussed in Section 3.4.

The concept of feature significance is central to the floating search methods. In this study, we extend this concept in the context of system identification and integrate the term orthogonalization to determine the significance of a term/term subset. This is discussed in Section 3.1. Further, each search method and associated adaption are discussed in the context of the system identification in Section 3.2, Section 3.3 and Section 3.4. Note that the adapted versions of SFFS, IFFS and OS are referred to respectively as Orthogonal Sequential Floating search (OSF), Orthogonal Improved sequential Floating search (OIF) and Orthogonal Oscillating Search (O2S).

3.1 Term Significance

Determination of feature significance is the core of the floating search methods [22, 23, 24, 25], where in each search step, the most significant feature/features are selected from the pool of available features and then the least significant feature/features are removed, from the selected feature subset. If a similar approach is to be applied to address the structure selection problem, it is essential to determine first the ‘significance’ of a term/term subset. For this purpose, in this study, we extend the definitions of feature significance developed by Pudil, Somol and co-workers [22, 23, 24] as follows:

Let $\mathcal{X}_{k}=\{x_{1},x_{2},\dots x_{k}\}$ be the term subset selected by the search algorithm at particular search step from the superset of NARX terms, $\mathcal{X}_{model}=\{x_{1},x_{2},\dots x_{n}\}$. Consequently, the subset of available terms is given by $\mathcal{X}_{model}\setminus\mathcal{X}_{k}$. Let $\mathcal{X}_{o}^{i}$ denote $i^{th}$ subset containing ‘$o$’ number of terms. The criterion function corresponding to the term subset $\mathcal{X}_{k}$ is denoted by $J(\mathcal{X}_{k})$ and is given by

[TABLE]

where, $ERR(x_{i})$ is given by (4).

For these subsets, the significance of term/term subset is defined as follows:

Definition 1**.**

The most significant term with respect to $\mathcal{X}_{k}$ is denoted as ‘$x^{MS}$’ and is defined as:

[TABLE]

Definition 2**.**

The least significant term in $\mathcal{X}_{k}$ is denoted as ‘$x^{LS}$’ and is defined as:

[TABLE]

Definition 3**.**

The most significant subset of ‘$o$’ number of terms with respect to $\mathcal{X}_{k}$ is denoted as ‘$\mathcal{X}_{o}^{MS}$’ and is defined as:

[TABLE]

Definition 4**.**

The least significant subset of ‘$o$’ number of terms in $\mathcal{X}_{k}$ is denoted as ‘$\mathcal{X}_{o}^{LS}$’ and is defined as:

[TABLE]

3.2 Orthogonal Sequential Floating Search

The principle of the generalized floating search was introduced by Pudil et al. in [22] which is essentially a sequential approach with a distinct ability to ‘backtrack’, i.e., the selected/discarded term can be discarded/selected in later stages. It is worth to note that the floating search can proceed either in forward (bottom-up) or backward (top-down) direction. In this study, we focus on the forward or bottom-up approach. This is referred to as Orthogonal Sequential Floating (OSF) search. The pseudo code for OSF is shown in Algorithm 1. Note that the OSF does not include ‘term swapping’ procedure (Line 1-1, Algorithm 1), in contrast to the improved floating search approach (which will be discussed in Section 3.3).

The search process in OSF begins with an empty term subset and is continued till the term subset with the desired cardinality ‘$\xi$’ is obtained. Each search step essentially involves two procedures: inclusion of the most significant term followed by the removal of the least significant term/terms through backtracking. This procedure is discussed briefly in the following:

Let $\mathcal{X}_{k}$ denote a term subset with $k$ number of terms at a particular search step. First, following Definition 1, the most significant term ($x^{MS}$) with respect to $\mathcal{X}_{k}$, is identified from the pool of available terms. This term is included in the term subset provided it leads to a better criterion function, i.e., $J(\mathcal{X}_{k}\cup x^{MS})>J(\mathcal{X}_{k+1})$. This procedure is outlined in Lines 1-1, Algorithm 1. After including a term, the focus is on removing least-significant terms present in the selected term subset through adaptive backtracking, as shown in Lines 1-1, Algorithm 1. It is worth to emphasize that the backtracking continues to remove the least significant term until an improvement in the criterion function is obtained, i.e., $J(\mathcal{X}_{k}\setminus x^{LS})>J(\mathcal{X}_{k-1})$.

To gain an insight into the search dynamics of the proposed OSF, consider the following system.

[TABLE]

For identification purposes, $1000$ input-output data-points are generated by exciting the system $\mathcal{S}_{1}$ with zero-mean uniform white noise sequence with unit variance, $u\sim WUN(0,1)$ and Gaussian white noise, $e\sim WGN(0,0.05)$. The model set of $165$ NARX terms is obtained with the following specifications: $[n_{u},n_{y},n_{l}]=[4,4,3]$.

The OSF is applied to identify six significant terms of $\mathcal{S}_{1}$, i.e., $\xi=6$. Note that the floating search methods require a priori specification of subset size, ‘$\xi$’. A simple approach to select the appropriate subset size $\xi$ and the associated issues will be discussed in Section 3.5. The search behavior of OSF is explained as follows: In each search step, the term subset, its cardinality and the corresponding criterion function are recorded, which are shown in Table 1 and Fig. 1. For the sake of clarity, the system terms and spurious terms are respectively denoted by ‘$x$’ and ‘$\mathcal{T}$’. Further, a better subset for the given cardinality is annotated by ‘+’, e.g., $\mathcal{X}_{2}^{+}$ indicates that $J(\mathcal{X}_{2}^{+})>J(\mathcal{X}_{2})$.

The positive effects of backtracking are clearly evident in Fig. 1. The backtracking occurs after the subset $\mathcal{X}_{5}$ is obtained. Note that the subsets obtained after backtracking ($\mathcal{X}_{2}^{+},\mathcal{X}_{3}^{+},\mathcal{X}_{4}^{+}$ and $\mathcal{X}_{5}^{+}$) yield improved criterion function in comparison to the previous subsets with the similar cardinality ($\mathcal{X}_{2},\mathcal{X}_{3},\mathcal{X}_{4}$ and $\mathcal{X}_{5}$). A closer examination of subsets given in Table 1 reveals that the backtracking could remove the spurious term ‘$\mathcal{T}_{1}$’ which is first included at Step 2 and remains in the selected subsets till Step 5. Further, the subset containing all the system terms, $\mathcal{X}_{5}^{+}=\{x_{1},x_{2},x_{3},x_{4},x_{5}\}$, is obtained during backtracking at Step 6. This subset is eventually selected at Step 9 when the forward inclusion procedure failed to yield a better subset, i.e., the search is restored to previous best subset (in this case $\mathcal{X}_{5}^{+}$, as $J(\mathcal{X}_{4}^{+}\cup x^{MS})<J(\mathcal{X}_{5}^{+})$; see Line 1-1, Algorithm 1).

3.3 Orthogonal Improved Floating Search

The next algorithm considered in this study is based on the improved floating search proposed in [25]. This algorithm adds a new procedure referred to as ‘swapping’ to replace a weak feature, besides retaining both the procedures of the floating search. The improved floating search is adapted for the structure selection by including the ERR metric associated with orthogonal least square algorithms of [4]. This is referred to as Orthogonal Improved Floating (OIF) search and discussed briefly in the following.

The procedures involved in OIF are shown in Algorithm 1. As discussed earlier, OIF retains the ‘forward inclusion’ (Line 1-1, Algorithm 1) and ‘backtracking’ (Line 1-1, Algorithm 1) procedures of the OSF and the same discussions are valid here. In addition, the swapping procedure is introduced after backtracking as shown in Line 1-1. The objective of this procedure is to replace a weak/non-significant term which is explained as follows:

Assume that $\mathcal{X}_{k}=\{x_{1},x_{2},\dots x_{k}\}$ is the term subset obtained after backtracking. To replace a non-significant term, several new candidate term subsets are generated from $\mathcal{X}_{k}$. The $i^{th}$ candidate term subset (denoted as $\hat{\mathcal{X}^{i}}$) is generated from $\mathcal{X}_{k}$ in the following two steps:

First, the $i^{th}$ term is removed from $\mathcal{X}_{k}$, i.e., $\hat{\mathcal{X}^{i}}\leftarrow\{\mathcal{X}_{k}\setminus x_{i}\}$ 2. 2.

Next, the most significant term with respect to $\hat{\mathcal{X}^{i}}$ is included, i.e., $\hat{\mathcal{X}^{i}}\leftarrow\{\hat{\mathcal{X}^{i}}\cup x^{MS}\}$

Following these steps, a total of ‘$k$’ candidate term subsets are obtained from $\mathcal{X}_{k}$ and the subset with the highest criterion function, $J(\cdotp)$, is referred to as ‘swapping subset’(denoted by $\hat{\mathcal{X}}$). If the swapping subset yields an improved criterion function, then it replaces the current term subset $\mathcal{X}_{k}$ and sent for backtracking, as outlined in Line 1-1, Algorithm 1.

Thus, OIF takes a two-pronged approach to remove weak/non-significant terms: the backtracking removes the weak term/terms and the swapping replaces a weak term. Especially, the swapping procedure encourages the exploration of new term subsets which is likely to yield improved search performance in comparison to OSF. To investigate this further, the OIF is applied to identify the significant terms of the system $\mathcal{S}_{1}$ considered in (10) under the similar test conditions.

The final subsets obtained in each step and the corresponding criterion function are shown in Table 2. The variations in the criterion function are shown in Fig. 2. Though, the overall search dynamics of both OSF (Fig. 1) and OIF (Fig. 2) may appear similar in nature, it is interesting to note that the backtracking in OIF occurs one step earlier (Step 5, Table 2) in comparison to OSF (Step 6, Table 1). This could be explained by the swapping procedure of OIF which enables the inclusion of missing system term ‘$x_{2}$’ one step earlier.

3.4 Orthogonal Oscillating Search

The third algorithm considered in this study is the oscillating search introduced by Somol and Pudil in [24] which is the natural successor of the floating search principle. It retains one of the key properties of floating search,i.e., the ability to backtrack. Unlike the floating search principle, where the focus is on the individual feature, the oscillating search focuses on a subset of features. The basic search engine in this method is the ‘swing’ procedure where the subset under consideration is perturbed by successive addition/removal and removal/addition of multiple features. In this study, the oscillating search is also adapted for the structure selection by the introduction of ‘orthogonalization’ and referred to as Orthogonal Oscillating Search (O2S). The steps involved in O2S are outlined in Algorithm 2 and discussed briefly in the following:

Let ‘$\xi$’ denote the required number of terms. The search is initialized with a subset, $\mathcal{X}_{\xi}$, containing $\xi$ number of terms. This initial subset is obtained by successively adding the most significant terms (See Definition 1) as outlined Line 2-2, Algorithm 2. Subsequently, this subset is perturbed throughout the search process by the ‘swing’ procedures. The rationale behind the swing procedure is that a better subset could be obtained by replacing multiple (say, ‘$o$’) number of weak or non-significant terms present in the current subset with the similar number of relevant or significant terms. This could be achieved by the following two ‘swing’ procedures: 1) Down Swing and 2) Up Swing.

During the ‘down swing’, at first ‘$o$’ number of least significant terms in $\mathcal{X}_{\xi}$ are removed which yields a subset of $(\xi-o)$ terms, $\mathcal{X}_{\xi-o}$. Subsequently, ‘$o$’ number of most significant terms are added to $\mathcal{X}_{\xi-o}$ to yield a new, and possibly a better subset of $\xi$ terms, $\hat{\mathcal{X}_{\xi}}$. This new subset is retained provided the criterion function is improved. The down swing procedure is outlined in Line 2-2, Algorithm 2.

The Up-swing procedure, as the name suggests, first adds ‘$o$’ number of most significant terms to $\mathcal{X}_{\xi}$ to yield a subset $\mathcal{X}_{\xi+o}$. Next, a new subset $\hat{\mathcal{X}_{\xi}}$ is obtained by removing the $o$ number of least significant terms from $\mathcal{X}_{\xi+o}$. The other aspects are similar to the down swing. This procedure is outlined in Line 2-2, Algorithm 2.

It is clear that both swing procedures require identification of the following two subsets with respect to the term subset under consideration: the subset of ‘$o$’ most significant terms (denoted as $\mathcal{X}_{o}^{MS}$; see Definition 3) and the subset containing ‘$o$’ weak/non-significant terms (denoted as $\mathcal{X}_{o}^{LS}$; see Definition 4). In this study, $\mathcal{X}_{o}^{MS}$ is identified using OSF search procedure described in Section 3.2. The set of weak terms $\mathcal{X}_{o}^{LS}$ is identified using the top-down or backward search variant of OSF.

Further, the ‘search depth’ (denoted by ‘$o$’) determines the number of terms which are either to be added or to be removed in a particular swing. The search depth is adaptive and its value ‘$o$’ is dependent on the search dynamics. The search begins with $o=1$ and if two successive swings fail to identify a better term subset then it is incremented by $1$. The search depth is reset to $1$ whenever swing leads to an improved subset. The search is terminated when ‘$o$’ reaches to a pre-specified maximum search depth, $\mathcal{O}_{max}$. Consequently, the depth of search can easily be controlled by varying $\mathcal{O}_{max}$. This is one of the distinct and important ability of O2S which regulates the search efforts as per the prevailing requirement.

To analyze the search behavior of O2S, it is applied to identify system terms of the system $\mathcal{S}_{1}$ considered in the numerical example (10) with the following specifications: $\xi=5$ and $\mathcal{O}_{max}=2$. During this search for system terms, the term subset and the corresponding criterion function obtained in each step are recorded and shown here in Table 3. The change in subset cardinality during the search process is shown in Fig. 3.

The search is initialized by forward inclusion of $5$ terms, as per Line 2-2, Algorithm 2. This initial subset selects one spurious term, instead of the system term $x_{2}$, as shown in Step 1, Table 3. The search begins by the down swing procedure which does not improve the subset, as seen in Step 3, Table 3. However, in the subsequent up-swing, a better term subset, $\mathcal{X}_{5}^{+}$, is identified, i.e., $J(\mathcal{X}_{5}^{+})>J(\mathcal{X}_{5})$. Note that at this stage, all the system terms have been identified and the spurious term is discarded, as seen in Step 5, Table 3. The subsequent down swing (Step 7) and up swing (Step 9) identifies the same subset with no improvement in the criterion function, $J(\cdotp)$. Consequently, the search depth is increased by ‘$1$’, i.e., $o=2$. The search terminates at Step 11 as the down swing could not find a better subset and any further increase in search depth will lead to $o>\mathcal{O}_{max}$.

Note that the variation in subset cardinality during the aforementioned swing procedures are clearly visible in Fig. 3. Especially, see the drop in subset cardinality at Step 10. The increase in search depth at this step requires the search for the term subset with $\{\xi-o\}=3$ number of terms.

3.5 Selection of Subset Cardinality (Model Order Selection)

The floating search algorithms require the specification of ‘subset cardinality’ or number of terms to be identified. Given that this is not known a priori, in this study, the following procedure is followed to estimate the subset cardinality or ‘model order’.

For the system under consideration, a search algorithm is applied to identify several term subsets of increasing cardinality in a predefined search interval, denoted by $[\xi_{min},\xi_{max}]$. A set of subsets thus obtained can be used to estimate the model order using an appropriate Information Criterion (IC). The objective here is to locate a ‘plateau’ or ‘knee-point’ in the information criteria which would indicate an acceptable compromise for the bias-variance dilemma. The term subsets corresponding to such knee-point or plateau can be selected as the system model.

To further understand this procedure, consider the system $\mathcal{S}_{1}$ in numerical example (10). For this system, OIF is applied to identify term subsets of increasing cardinality in the range of $[\xi_{min},\xi_{max}]=[2,20]$. Thus, a family identified term subsets, denoted by ‘$\Omega$’, is obtained as follows:

[TABLE]

Next, for each identified subset, $\mathcal{X}_{\xi}\in\Omega$, the various information criteria given in A are determined. Fig. 4 shows the variation in information criteria as cardinality is varied from $\xi_{min}$ to $\xi_{max}$. It is observed that the ‘knee-point’ for all the criteria is obtained at $\xi=5$. This coincides with the actual cardinality of the system $\mathcal{S}_{1}$.

In several comparative investigations [26, 18], the Bayesian Information Criterion (BIC) was found to be comparatively robust among the existing information criteria. Therefore, for the remainder of the study, the cardinality corresponding to the minimum BIC is selected as the model order. Although following this approach the correct model order could be obtained for all test system (as will be shown in Section 5.1), in practice, we recommend the closer inspection of several term subsets surrounding the plateau or knee-point; as the existing research indicated that the information criteria tend to over-fit [26].

Since the existing research suggest that the dynamics of many nonlinear systems can be captured with $15$ or fewer significant terms [10], a pragmatic approach is followed in this study to fix the search interval of the cardinality as follows: $[\xi_{min},\xi_{max}]$=$[2,20]$. This can be considered as a thumb-rule. However, some exception may exist in certain non-linear systems. A simple procedure to determine the cardinality interval in such a scenario is discussed in Section 5.2.

4 Investigation Framework

In the following, the framework of this investigation is discussed. The overall procedure followed for the structure selection is shown in Algorithm 3. The test nonlinear system considered in this study are discussed in Section 4.1. Further, the possible search outcomes are discussed in Section 4.2, to evaluate the term subset found by the search algorithms qualitatively.

4.1 Test Non-linear Systems

In this study, a total of $6$ non-linear system have been selected from the existing investigations on structure selection [27, 12, 28, 13, 6, 7] and shown here in Table 4. The systems are excited by a white noise sequence having either uniform or Gaussian distribution as shown in Table 4. For identification purposes, a total of $1000$ input-output data points, $(u,y)$, are generated from each system. The structure selection is performed following the principle of cross-validation; where $700$ data points are selected for the estimation purpose and the remaining data points are used for validation, i.e., $\mathcal{N}_{v}=300$. For each system, the NARX model is generated by setting the input-output lags and the degree of non-linearity to: $[n_{u},n_{y},n_{l}]=[4,4,3]$. This gives the NARX model set, $X_{model}$, with a total of $165$ terms following (2), i.e., $n=165$.

4.2 Search Outcomes

For comparative evaluation purposes, the term subset found by each search algorithm is qualitatively evaluated. For this purpose, the following term sets are defined with reference to the NARX model given by (3),

• •

$\mathcal{X}_{model}$ : the set containing all terms of the NARX model

• •

$\mathcal{X}^{\star}$ : the optimum term subset or the set of system terms, $\mathcal{X}^{\star}\subset\mathcal{X}_{model}$

• •

$\mathcal{X}^{alg}$ : subset of terms identified by the search algorithm, $\mathcal{X}^{alg}\subset\mathcal{X}_{model}$

• •

$\mathcal{X}_{spur}$ : set of spurious terms which are selected by the search algorithm, but are not present in the actual system, i.e., $\mathcal{X}_{spur}=\mathcal{X}^{alg}\setminus\mathcal{X}^{\star}$

• •

$\varnothing$ : the null set

On the basis of these definitions, one of the following four search outcome can be expected:

Identification of the Correct Structure (Exact Fitting) :

In this scenario the identified model contains all the system terms and does not include any spurious terms, i.e., $\mathcal{X}^{alg}=\mathcal{X}^{\star}$ and $\mathcal{X}_{spur}=\varnothing$ 2. 2.

Over Fitting :

The identified model contains all the system terms; however spurious terms are also selected, i.e., $\mathcal{X}^{alg}\supset\mathcal{X}^{\star}$ and $\mathcal{X}_{spur}\neq\varnothing$ 3. 3.

Under Fitting-1 :

The algorithm fails to identify all the system terms; though it does not include any spurious terms, i.e., $\mathcal{X}^{alg}\subset\mathcal{X}^{\star}$ and $\mathcal{X}_{spur}=\varnothing$ 4. 4.

Under Fitting-2 :

The algorithm fails to identify all the system terms; however spurious terms are selected, i.e., $\mathcal{X}^{alg}\not\supset\mathcal{X}^{\star}$ and $\mathcal{X}_{spur}\neq\varnothing$

Thus, qualitatively, the search is successful when all the system/significant terms are identified, i.e., when the outcome is either Exact-Fitting or Over-Fitting. Note that at this stage, the structure identified by the algorithms are unaltered. The inclusion of few spurious terms can be tolerated provided that all significant terms are included (i.e. Over-Fitting scenario) as the spurious terms can easily be identified and removed through a simple null-hypothesis test on the corresponding coefficients.

5 Results

The goal of this study is to investigate the suitability of existing feature selection algorithms for the task of non-linear system identification. For this purpose, 3-well known floating search algorithms have been adapted in the proposed orthogonal floating search framework: OSF, OIF and O2S, which are discussed in Section 3. The adapted algorithms are initially applied to identify the significant terms of $6$ test non-linear systems described in Section 4.1. The search performance of the algorithms is compared on these systems from various perspectives in Section 5.1. The issues related to selection of cardinality interval $[\xi_{min},\xi_{max}]$ are discussed in Section 5.2.

Further, it is worth to note that while the proposed orthogonal floating search can work with any suitable criterion function, in this study, the Error-Reduction-Ratio (ERR) has been selected for this purpose due to its relative simplicity. However, ERR is often criticized for its ‘local’ nature and blamed for the unsatisfactory search outcome. This issue is investigated via a numerical example in Section 5.3.

Finally, the proposed orthogonal floating algorithms are applied to identify a discrete model of a continuous time system. The identified discrete model is validated through generalized frequency response functions. This part of the investigation is discussed in Section 5.4.

5.1 Comparative Evaluation

The objective of this part of the study is to compare the performance of orthogonal floating search algorithms. Following the procedure outlined Algorithm 3, the significant terms of the system shown in Table 4 are identified. It is worth to note that while OSF and OIF do not have any control parameter, O2S requires the selection of maximum search depth, $\mathcal{O}_{max}$. In this study, a maximum of $25\%$ terms is allowed to be exchanged, i.e., $\mathcal{O}_{max}=\lceil 0.25\times\min\{\xi,n-\xi\}\rceil$.

The outcomes of model order selection are shown in Fig 5. Note that the ‘knee-point’ in BIC coincides with the known cardinality (actual number of terms) of the systems. The selected cardinality (number of terms) for each system is shown in Table 5 along with the qualitative search outcomes. It is interesting to see that the best possible search outcome, ‘Exact-Fitting’, is obtained for all combinations of system and search algorithm, except one. For the system $\mathcal{S}_{1}$, an over-fitted model (with one spurious term) is selected by OSF as seen in Fig. 5(a) and Table 5. Note that this spurious term can easily be removed via a simple null-hypothesis test. Nevertheless, it is pertinent to identify the cause of this ‘over-fitting’ which could be ascribed either to the information criterion or to the search algorithm.

To investigate this further, the term selection behavior of the algorithms is determined. In particular, as the cardinality, $\xi$, is increased progressively from $\xi_{min}$ to $\xi_{max}$, the selection frequency of system terms is observed over the subsets identified by the algorithm, as follows:

[TABLE]

where, ‘$\tau_{i,\xi}$’ encodes whether the $i^{th}$ term, ‘$x_{i}$’, is included in the term subset $\mathcal{X}_{\xi}$ identified by the algorithm. Following this procedure, the term selection behavior is determined for all algorithms and shown here in Fig. 6, 7 and 8.

We first focus on the ‘over-fitting’ outcome by OSF on $\mathcal{S}_{1}$. The term selection by OSF for $\mathcal{S}_{1}$ is depicted in Fig. 6(a). It is observed that all the system terms are first identified at $\xi=6$ which is slightly higher than $\xi^{\star}$. This explains the ‘over-fitting’ outcome. For the same system, OIF (Fig. 6(b)) and O2S (Fig. 6(c)) could identify all system terms from $\xi=\xi^{\star}$. This further underlines the efficacy of the ‘term swapping’ procedure introduced in OIF.

Further, although the BIC could identify correct model order for all the systems in this study, the possibility of ‘over-fitting’ cannot be ruled out. We therefore focus on the search behavior for $\xi>\xi^{\star}$. It is easier to follow that all the system terms must be included in the term subsets with cardinality $\xi\geq\xi^{\star}$. Thus it is desirable to have $\tau_{i,\xi}=1,\forall\ \xi\geq\xi^{\star}$. If this condition is satisfied then the correct system structure can still be obtained from the over-fitted structures through an appropriate secondary refinement procedure, e.g., null-hypothesis test on the coefficients. The results show that this condition is satisfied for all systems by OIF and O2S, as seen in Fig. 6, 7 and 8. OSF could also satisfy this condition for all the systems except $\mathcal{S}_{1}$.

5.2 Comments on the Interval Selection

The objective of this part of the study is to discuss the issues related to the selection of the cardinality interval. Given that the cardinality of the system under consideration is not known a priori, a careful selection of the cardinality interval is crucial. In the previous subsections, a pragmatic approach is followed to select the cardinality interval, i.e., $[\xi_{min},\xi_{max}]=[2,20]$; since in practice, over-parameterized models with higher cardinality are usually not desirable.

However, the selection of cardinality interval is system dependent. Therefore, the cardinality interval should be selected judiciously, especially for the higher order systems. In order to gain further insight into this issue, consider the following nonlinear system with $23$ terms,

[TABLE]

For the identification purposes, $1000$ input-output data points are gathered with $u\sim WUN(0,1)$ and $e\sim WGN(0,0.01^{2})$. A total of $286$ candidate terms are obtained with the following specification for the NARX model: $[n_{u},n_{y},n_{l}]=[5,5,3]$.

Initially, the system model is identified following the similar procedure, where the cardinality interval is specified as: $[\xi_{min},\xi_{max}]=[2,20]$. The outcome of the cardinality selection for this interval is shown in Fig. 9(a) which shows that the BIC is minimum at $\xi=20$. However, it is not certain whether this is the desired ‘knee-point’. Therefore, the cardinality interval is increased to $\xi_{max}=50$. The consequent outcome is shown in Fig. 9(b). This result indicates a continual increase in BIC for $\xi\geq 30$ and thereby confirms that the plateau around $\xi\in[22,29]$ does contain the desired knee-point.

As illustrated by this numerical example, the cardinality interval can easily be adjusted by observing the consequent ‘$\xi$ vs. BIC’ relationship. We recommend the OSF to determine/adjust the cardinality interval due to its inexpensive computational requirements.

5.3 Orthogonal Floating Search Under Slow Varying Input

As discussed in Section 3.1, the present investigation uses the Error-Reduction-Ratio (ERR), instead of the other performance metrics, such as simulation error, ultra least square criterion. However, in recent years, the OFR-ERR has been criticized for leading to suboptimal search performance under certain scenarios [13, 15, 7, 6], where it is argued that this is due to the ERR criterion.

However, the results of this investigation on several benchmark nonlinear systems convincingly demonstrate that ERR could be an effective metric, provided it is paired with an appropriate search strategy. Thus, it can be argued that the suboptimal performance of OFR-ERR may be the consequence of the search strategy and not of the performance metric. Note that OFR-ERR is essentially a unidirectional sequential search approach which leads to the ‘nesting effect’, i.e., once a term is included it cannot be discarded. Thus the search performance may significantly be improved if the nesting effect could be alleviated. The proposed orthogonal search framework effectively addresses this issue and therefore it can yield better search performance even with a relatively simple metric such as ERR.

To demonstrate this, we consider the following numerical example, which has been used by various researchers for comparative investigation [13, 15, 16]:

[TABLE]

where, ‘$y(\cdotp)$’ is the observed value of ‘$w(\cdotp)$’. The system is excited with the slow-varying input $u(\cdotp)$ which is an AR process with a repeat pole at 0.8. A total of 165 candidate terms are obtained by the following NARX specifications: $[n_{u},n_{y},n_{l}]=[4,4,3]$.

The system is identified by the following three algorithms: OFR-ERR, OIF and O2S. The search dynamics of the compared algorithms are recorded and are shown in Table 6 (OFR-ERR), Table 7 (OIF) and Table 8 (O2S). Note that, initially, all the algorithms selected the autoregressive terms ‘$y(k-1)$’ and ‘$y(k-2)$’. As seen in Table 6, OFR-ERR could not remove the autoregressive terms due to the nesting effect. Furthermore, several other spurious terms have been selected instead of two system terms: ‘$u(k-2)$’ and ‘$u(k-1)u(k-2)$’. In contrast, OIF could easily remove the autoregressive terms through swapping procedure (see Step 3 and 4, Table 7). Similarly, the autoregressive terms are replaced by the system terms and the correct structure is identified in the first down-swing procedure of O2S (see Step 3, Table 8).

It is obvious that the floating framework enabled OIF and O2S to identify a better subset for the given cardinality. For example, all the system terms are identified by OIF and O2S at $\xi=4$ with the criterion function $J(\cdotp)=0.99757$ (see Table 7-8). In contrast, OFR-ERR identified only one system term at $\xi=4$. This is further reflected in the relatively lower value of the criterion function, $J(\cdotp)=0.98738$ (see Table 6). This example thus clearly demonstrate that ERR could be an effective performance metric when it is paired with a proper search strategy.

It is worth to emphasize that, for identification purposes, the input signal to the system should be persistently exciting in nature. The example considered in (14) is the exception to this fundamental fact, and it represents a worst-case scenario. In practice, such data should not be used for identification purposes. This example has been considered only for comparative analysis.

5.4 Discrete models of Continuous-time Systems

The effectiveness of the proposed framework is further demonstrated by fitting a polynomial NARX model to a continuous-time system. Since a discrete model of a continuous time system is not unique, the identified discrete-time model is validated by comparing Generalized Frequency Response Functions (GFRFs). Note that if the identified model has captured the system dynamics, the corresponding GFRFs should match.

For this purpose, the Duffing’s Oscillator is considered, which is given by,

[TABLE]

where, $\omega_{n}=45\pi$, $\zeta=0.01$ and $\varepsilon=3$. The input-output data are generated by exciting the system with $u(\cdotp)\sim WUN(0,1)$ and the output is sampled at $500\ Hz$. A total of $1000$ data points are generated. The discrete NARX model is fitted using $700$ data points, and the remaining 300 data points are used for validation. The super-set containing a total of $286$ NARX terms is obtained by the following specifications: $[n_{u},\ n_{y},\ n_{l}]=[5,\ 5,\ 3]$.

The following model is identified by all the proposed algorithms OIS, OIF and O2S:

[TABLE]

The fitting error obtained over the validation data points is $3.96\times 10^{-6}$.

The linear and the third order GFRFs of the continuous system are shown in Fig. 10(a), 10(c) and 10(e). As indicated by these results, the frequency response of the Duffing’s oscillator is characterized by the presence of several ‘ridges’ and ‘peaks’. The ridges occur when the system is excited by the linear resonant frequency (i.e., $f_{1}$, $f_{2}$ or $f_{3}=\pm 22.5$ Hz) or when the following condition is satisfied: $f_{1}+f_{2}+f_{3}=\pm 22.5$. These ridges are shown by ‘dotted’ lines in the contour plot shown in Fig. 10(e). Further, the peaks occur whenever all frequencies are equal to $\pm 22.5$. Note that these observations are in line with the earlier investigation [29].

The GFRFs obtained with the discrete identified model are shown in Fig. 10(b), 10(d) and 10(f). It is clear that the identified model could capture all the characteristics features of the Duffing’s oscillator.

6 Discussion & Conclusions

A new Orthogonal Floating Search Framework has been proposed to identify significant terms of polynomial NARX models. The proposed framework integrates the classical ERR metric with floating search algorithms which retains the simplicity of OFR while eliminating its nesting effects. Following this framework, three well-known feature selection algorithms, such as SFFS, IFFS and OS, have been adapted to solve the structure selection problem. The effectiveness of this approach has been demonstrated considering several benchmark nonlinear systems including a system excited by a slowly varying input. The later system is often considered as a worst-case scenario. Further, a discrete-time model is fitted to Duffing’s oscillator and it is validated by comparing generalized frequency response functions. These results convincingly demonstrate that the floating feature selection algorithms can easily be tailored for the purpose of structure selection in nonlinear system identification. For example, relatively simple procedures of the floating search, such as backtracking, term swapping and swing, can effectively counter the nesting effect and thereby yield significant improvement in the search performance (as shown in the results). The results of this investigation support the argument that the suboptimal performance of OFR-ERR is due in major part to its poor search strategy.

Appendix A Information Criteria

Let $\mathcal{X}_{\xi}$, $\xi\in[\xi_{min},\xi_{max}]$ denote a term subset identified the search algorithm. For this subset and the corresponding model-predicted output ($\hat{y}$), the information criteria are given by,

• •

Akaike Information Criterion (AIC):

[TABLE]

• •

Bayesian Information Criterion (BIC):

[TABLE]

• •

Final Prediction Error (FPE):

[TABLE]

• •

Khundrin’s Law of Iterated Logarithm Criterion (LILC):

[TABLE]

where,‘$\mathcal{N}_{v}$’ denotes the length of the validation data; $\mathcal{E}=\frac{1}{\mathcal{N}_{v}}\sum\limits_{k=1}^{\mathcal{N}_{v}}[y_{k}-\hat{y}_{k}]^{2}$