Two-Dimensional (2D) Particle Swarms for Structure Selection of Nonlinear Systems

The present study proposes a new structure selection approach for non-linear system identification based on Two-Dimensional particle swarms (2D-UPSO). The 2D learning framework essentially extends the learning dimension of the conventional particle swarms and explicitly incorporates the information about the cardinality, i.e., number of terms, into the search process. This property of the 2D-UPSO has been exploited to determine the correct structure of the non-linear systems. The efficacy of the proposed approach is demonstrated by considering several simulated benchmark nonlinear systems in discrete and in continuous domain. In addition, the proposed approach is applied to identify a parsimonious structure from practical non-linear wave-force data. The results of the comparative investigation with Genetic Algorithm (GA), Binary Particle Swarm Optimization (BPSO) and the classical Orthogonal Forward Regression (OFR) methods illustrate that the proposed 2D-UPSO could successfully detect the correct structure of the non-linear systems.


Introduction
Construction of mathematical models from observed input-output data; popularly known as system identification, has been a major research concern from diverse fields such as statistics, control theory, information theory, economics, ecology, and agriculture. Most of the practical systems are inherently nonlinear and therefore the development of system identification methods, which can be valid for a broad class of nonlinear systems, has attracted many researchers [1][2][3][4].
In particular, this study aims to develop a new approach for the identification of nonlinear systems represented by polynomial Nonlinear Auto-Regressive with eXogenous inputs (NARX) models [5][6][7]. The NARX model enables the representation of a broad class of nonlinear systems in recursive lagged input and output terms which leads to a convenient linear-in-parameter form. The major challenge of the non-linear system identification is to determine which input-output terms from the NARX model shall be included in the system model to capture its dynamics, referred to as the structure selection problem. This difficulty can be ascribed to the exponential increase in the number of terms with the increase in the degree of non-linearity and maximum lags of inputs and outputs of the NARX model.
In the sequential search, a structure is built incrementally by adding and/or removing a term in each step. Among the sequential search methods, Orthogonal Forward Regression (OFR) with Error Reduction Ratio (ERR) has been extensively studied and applied to many applications by Billings and co-workers [16,[23][24][25][26]. OFR is essentially a sequential greedy search method, wherein the structure is built incrementally by including the term with the highest performance metric (i.e., ERR) in each step. The main limitation of OFR-ERR is the emphasis on a single term; though the terms are decoupled through orthogonalization, the performance metric such as ERR may depend on the order in which terms are orthogonalized [27]. Several variants have been proposed to further improve the performance of OFR [17,18,[27][28][29]. Recently, an alternate approach referred as 'RJMCMC', is proposed [11] in which a term is either added through 'birth move' or removed through 'death move'. The probabilities for these moves are determined using the Bayesian inference.
In contrast, in the subset search, entire candidate structure/term subset is determined and evaluated in each step. The common element in all subset search methods is the use of a stochastic component to enable better sampling of search space. Earlier research in this direction includes the application of the evolutionary algorithms such as Genetic Algorithm (GA) and Genetic Programming (GP) [19,30,31]. In recent years, various other structure selection algorithms have been proposed which also fall into the subset search category, e.g., see [9,10,12,32]. In the randomized search proposed in [9], known as 'RaMSS', each term is assigned with a selection probability. A term is included to/excluded from the term subset on the basis of the corresponding probability. In each iteration, a fixed number of term subsets are generated and evaluated. On the basis of this evaluation, the selection probability of each term is updated.
It is worth to mention that the goal of the structure selection is to identify a model which can capture the linear and non-linear dynamics of the underlying system while being parsimonious enough to be practical [3,[33][34][35]. The major challenge for any search method is to balance these contradictory objectives, i.e., bias vs. variance trade-off. The search for optimum structure, therefore, entails selection of correct/significant terms as well as the number of terms.
This study proposes a new approach to model structure selection, where the information about the terms and the number of terms (referred here as cardinality) is explicitly integrated into the search algorithm to determine the correct structure. The core idea is to use the cardinality/ the number of terms as the additional information to balance bias-variance dilemma. The proposed approach is, in essence, a population based search heuristic which has been developed in particle swarm theory [36], where each search agent (particle) encodes a term subset/structure. In addition, each particle independently stores the selection likelihoods of both, the cardinality (number of terms) and the terms. In each iteration, the best structures, obtained hitherto, are used as learning exemplars. Further, a novel learning framework is developed to extract the information about the beneficial/significant terms and cardinality from the learning exemplars/best structures. This information is used to update the selection likelihoods of the terms and cardinality. Subsequently, a new particle/structure is determined through the inclusive use of both the selection likelihoods. In essence, the proposed approach extends the learning dimension of the conventional particle swarms to integrate information about cardinality. Due to this distinctive property, it is referred as 'Two-Dimensional (2D)' learning framework.
Note that the 2D-learning was originally developed by the authors for the 'feature selection' problem in machine learning [37]. The proof of concept of this algorithm as an alternate method of structure selection of nonlinear systems has been reported in [38]. The present investigation significantly differs from [38] in the following aspects: • The proposed algorithm is better formalized and explained through illustrative example from the perspective of system identification. • A detailed analysis of the search behavior and the influence of the control parameters on the search performance has been reported. • The comparative evaluation is carried out considering 4-established structure selection approaches. • A rigorous evaluation of the proposed approach has been carried out considering 7 benchmark nonlinear systems which are taken from the existing research. • The search behavior of the algorithm, for a system with slowly varying excitation, has been investigated. This may be considered as a probable worst case scenario. • A discrete time model of the continuous time system has been identified and validated through Generalized Frequency Response Functions (GFRF) [39]. • In addition to the simulated examples, a practical case study on identification of non-linear wave-forces has been presented. The rest of the article is organized as follows: the NARX model and the structure selection problem are briefly discussed in Section 2.1. The philosophy of the 2D learning is discussed in Section 3. The framework of this study is described in Section 4. The results are discussed at length in Section 5, followed by the conclusions in Section 6.

The Structure Selection Problem
The first step in the nonlinear system identification is the selection of system representation such as Volterra, NARX, Wiener, neural networks, polynomial models, and rational models. Among these, the NARX models yield the simplest system representation while providing precise information about system dynamics. The focus of this study is, therefore, the identification of nonlinear systems represented by polynomial nonlinear auto-regressive with exogenous inputs (NARX) models [5,6]. In the following, the NARX model and the structure selection problem are briefly discussed.

The Polynomial NARX Model
The NARX model represents a non-linear system as a function of recursive lagged input and output terms as follows: 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 {· } is some nonlinear function of degree n l .
The total number of possible terms or model size (N t ) of the NARX model is given by, This model is essentially linear-in-parameters and can be expressed as: where, x 1 (k) = 1, and p y , q u ≥ 0; 1 ≤ p y + q u ≤ n l ; 1 ≤ n y j ≤ n y ;1 ≤ n u k ≤ n u ; n l is the degree of polynomial expansion.

Problem Formulation
Consider the identification of a nonlinear system represented by polynomial NARX model. Given a large model set with N t number of terms, denoted as, where, x 1 , x 2 , . . . x Nt could represent any possible linear or non-linear term of the NARX model. The goal of the structure selection is to determine the correct/optimum subset of terms, X ⊂ X model , by minimizing a suitable criterion function,'J(· )', i.e., where, 'ξ X ' denotes the number of terms present in term subset 'X'. This will be referred as 'cardinality' in this study. Note that, this is essentially a combinatorial optimization problem. An exhaustive search of all possible term subsets to solve this problem is often intractable even for a moderate number of terms N t , as it requires the examination of 2 Nt term subsets. Since the information about cardinality (number of terms) is not known a priori, a successful structure selection method should address the following: 1) How many terms are required to represent the system dynamics? and 2) Which terms should be included? These two issues are crucial to effectively address the bias-variance dilemma.
However, most of the existing search methods focus on evaluating the significance of the terms. To the best of our knowledge, the information about the number of terms (i.e., cardinality) has not been exploited for the benefit of the search process. This study, therefore, proposes a new approach with the following key features: • The information about the subset cardinality (number of terms) is determined and explicitly integrated into the search process. Instead of assigning the selection likelihoods only to the terms, these are assigned both to the number of terms (cardinality) and the terms. These likelihoods are updated continually throughout the search process. The new structures are explored through the joint use of these selection likelihoods. • In most of the existing search approaches, a common selection probability is used to generate test structures. In contrast, in the proposed approach, each search agent/particle independently stores and updates the selection likelihoods. Thus, each particle can explore structures having different combinations of cardinality and the terms. This, in theory, often translates into the improved sampling of the search space.

Proposed Two-Dimensional (2D) Algorithm for Structure Selection
The Particle Swarm Optimization (PSO), which is the kernel of the 2D learning, is based on the concept of search through the social co-operation. In essence, PSO employs a group (swarm) of mass-less objects, referred to as 'particles', to represent a possible solution to the problem under investigation. Each particle has two memory attributes: 1) the best solution found by itself (cognitive memory) and 2) the best solution found by the entire swarm (social memory). These memory attributes are updated throughout the search process and used as 'learning exemplars', i.e., the memory attributes are used to determine a new move of each particle on the search landscape. Further details about the particle swarm optimization, its discretization framework and the 2D learning framework can be found in [36,37,40].
The proposed 2D learning framework has been developed based on the particle swarm theory, and this is applied to determine the correct structure of a nonlinear system represented by the polynomial NARX model. In this study, we focus on the following two key issues: 1) How many terms should be included in the structure (to be called as cardinality in the present context)? and 2) Which terms should be included ? For this purpose, the learning dimension of particles is extended to explicitly integrate the cardinality information into the search process, which results in the proposed 2D learning framework, as discussed in the following.

Philosophy of 2D Learning
To understand the Two-Dimensional (2D) learning in the context of system identification, consider the problem of selecting the correct structure from the NARX model with 'N t ' number of terms. Let the set containing all model terms be denoted as X model = x 1 x 2 . . . x Nt . For this problem, each particle (β) encodes a candidate structure or set of terms in an N t -dimensional binary vector. For example, the i th particle, β i , is represented as: The m th term (x m ) from X model is included into the candidate structure provided the corresponding bit in the particle, 'β i,m ' is set to '1'. For more details, see the illustrative example in Appendix A.
The key feature of the 2D learning is the explicit integration of information about the cardinality (the number of terms, denoted as 'ξ') into the search process. In this learning framework, the selection likelihoods of both the cardinality and terms are independently stored in the 'velocity' (V ) of each particle. For this purpose, the learning dimension of particles is extended, and the velocity of the i th particle is represented as a two-dimensional matrix of size (2 × N t ) as follows: Each row of 'V ' corresponds to a separate learning dimension. The elements in the first row of V store the selection likelihoods of cardinality (i.e., how many terms to be included ), and the elements in the second row give the selection likelihoods of the corresponding terms (i.e., which terms to be included ). For example, 'v i 1m ' gives the probability of including a total of 'm' number of terms in the i th structure, i.e., ξ i = m. On the other hand, 'v i 2m ' gives the probability of including the m th term in the i th structure.
It is worth to emphasize that the 'velocity' is the key search engine of the 2D learning. Throughout the search process, the cardinality and term selection likelihoods, stored in the velocity, are continually updated using the information extracted from the learning exemplars. In essence, each iteration of the 2D learning involves the following three steps: Re-initialize the velocity of the particle 10 Set count i to zero Evaluate the learning sets, L α1 , L α2 , L α3 and L i , as per Algorithm 2 13 Update the velocity of the i th particle as per (5) 14 Update the position of the i th particle following Algorithm-3 The ring-topology is used to define the particle neighborhood 1. Evaluation of the Learning Sets: Beneficial information about the cardinality and the terms is extracted from the learning exemplars and stored in the learning sets. 2. Velocity Update: The learning sets are used to update the cardinality and term selection likelihoods. 3. Position Update: A new structure/term subset is determined through joint use of the cardinality and term selection likelihoods. The following subsections describe these procedures in detail. The overall procedure involved in 2D-UPSO is described in Algorithm 1. The time complexity of this algorithm has been studied and reported in [37] where it has been shown that the 2D-UPSO could complete the search in comparatively shorter run time.

Evaluation of Learning Sets
The core idea of the particle swarm theory is to learn from the memory attributes/learning exemplars to determine better solutions. In 2D learning, this is accomplished through the extraction of learning sets from the learning exemplars. To understand this process, let 'α' denote a memory attribute/learning exemplar such as 'personal best' (pbest), 'swarm best' (gbest) or 'neighborhood best' (nbest). Note that, α is essentially the best structure (term subset) found hitherto in the search process. Thus, it is an N t -dimensional binary string.
The objective of the 2D learning is to extract the following information from the learning exemplar 'α': 1) cardinality and 2) the terms that have been included in 'α' but not in the particle, 'β'. Note that the information about the cardinality and the terms is stored independently into a two-dimensional binary learning set, L, of size (2 × N t ). A typical learning set is thus given by, Similar to V , the first row of L stores the learning about the cardinality and the second row stores the learning about beneficial terms acquired from the learning exemplar.
Given that both 'α' and 'β' are essentially binary strings, the required information to determine the learning set (L) can easily be extracted following the steps outlined in Algorithm 2. See Appendix B for an illustrative example. The similar procedure is followed to extract a learning set from each learning exemplar such 'pbest' and 'gbest'.

Velocity Update
Once the learning sets have been extracted, the velocity of each particle is updated through a simple velocity update rule. For example, at any given iteration t, the velocity of the i th particle is updated as: where u f ∈ [0, 1] denotes the unification factor ; r 1 , r 2 , r 3 ∈ R 2×Nt are uniformly distributed random numbers in [0, 2]. Further, L α 1 , L α 2 , L α 3 and L i denote respectively the learning set derived from pbest, gbest, nbest and β i . One of the key features of 2D learning is the inclusion of the self-learning set, L i , which is derived from the particle position, β t i . The influence of this learning set is controlled through '∆' which is defined as follows: where, 'J t i ' and ' − → J t ' denote respectively the criterion function corresponding to the i th particle and the entire swarm at iteration t. The objective here is to provide a positive feedback, 'δ i ', if the particle position, β t i , deemed beneficial. Further, the value of δ is dependent on the performance of a particle, e.g., the particle with the minimum fitness 'min( − → J t )' is assigned with the maximum δ, as seen in (6).
Note that, in the proposed 2D learning framework, the velocity update rule in (5) uses Unified Particle Swarm Optimization (UPSO) [41], which is a well-known PSO variant. Therefore, it is referred to as '2D-UPSO' throughout this manuscript.

Position Update
The position update procedure is carried out in two stages, to integrate the information about the cardinality and the terms, as outlined in the Algorithm 3. In the first stage, the cardinality of a new structure is determined. Subsequently, the required beneficial/significant terms are selected. This process is further explained through an example in Appendix C.

Criterion Function
As discussed earlier, each search agent in the population based meta-heuristics represents a particular term subset (e.g., a parent in GA, a particle in BPSO and 2D-UPSO). The criterion function, 'J(· )', is essentially a performance metric which estimates the quality of the search agent/term subset under consideration. For most of the structure selection problems, the error in model-predicted output (ŷ), over the validation data, is often considered to be a good performance metric [42]. Further, it is equally important to remove the spurious terms to prevent over-fitting of Algorithm 3: The position update of the i th particle in 2D-learning Input : v i ; velocity of the i th particle Output: β i ; position of the i th particle Selection of the cardinality based on roulette wheel, (ξ i ) the model. Hence, the structure selection problem is multi-objective in nature wherein the objective is to find a parsimonious term subset with the minimum fitting error. To balance these objectives, several information-theoretic criteria have been proposed such as Akaike Information Criterion, Bayesian Information Criterion (BIC), Minimum Description Length and the others [42,43]. A comparative analysis of the information criteria suggests that BIC is robust in the presence of various levels of measurement noise [42]. It is therefore selected as the criterion function in this study. For a given term subset, X i , this is defined as: where, 'E' denotes the error with respect to the model-predicted output (ŷ) which is obtained with X i ; 'N v ' denotes the length of the validation data and 'ξ X i ' is the cardinality (number of terms) in 'X i '.

Investigation Framework
The following subsections provide the details about the framework of this investigation. In this study, the search performance of 2D-UPSO is compared with the conventional Binary Genetic Al- Record the best term subset (X k run ) and the corresponding criterion function (J(X k run )) found by the algorithm  [44], Orthogonal Forward Regression (OFR-ERR) [24] and OFR with pruning [45]. For this purpose, the search environment is set as per Section 4.1. The search performance of all algorithms is evaluated through 6 test nonlinear systems which are described in Section 4.3. Further, the possible search outcomes are discussed in Section 4.2, to evaluate the term subset found by the search algorithms qualitatively.

Search Setup
The experimental setup of this investigation is described in Algorithm 4. To accommodate the stochastic nature of GA, BPSO and 2D-UPSO, 'R' number of independent runs of each algorithm are carried out on the test systems (R is set to 40). Each run is set to terminate after 6000 Function Evaluations (FEs) of criterion function, J(· ), given by (7). In each run, the best term subset and the corresponding J(· ) found by the algorithm are recorded. Due to the stochastic nature of the search algorithms, the best term subset from R independent runs (X alg ) is selected, as outlined in Algorithm 4, Line 4-9.
Note that a judicious selection of search parameters is crucial to the performance of any search algorithm. The proposed 2D-UPSO has two search parameters: unification factor (u f ) and refresh gap (RG). The selection of these parameters and their effects on the search performance have been discussed in Section 5.1. The search parameters of GA and BPSO have been selected on the basis of previous research [44,46] as follows: • Further, most of the existing structure selection methods are likely to select some terms which are not present in the model, (spurious terms). Although the number of spurious terms chosen by proposed 2D-UPSO is significantly less compared to GA and BPSO (as will be shown in Section 5.2), the selection of additional terms cannot always be ruled out. A simple null-hypothesis test may be applied to the coefficients to remove such spurious terms and refine the identified structure further.

Search Outcomes
For comparative evaluation purposes, the term subset found in each run of the compared algorithm is qualitatively evaluated. Given that the proposed 2D-UPSO and the compared algorithms (GA and BPSO) are stochastic in nature, the objective of the qualitative comparison is to evaluate the robustness of the search performance.
For this purpose, the following term sets are defined in reference to the NARX model given by (2), • X model : the set containing all terms of the NARX model • X : the optimum term subset or the set of system terms, X ⊂ X model • X k run : term subset found in the k th run of the algorithm, k = 1, 2, . . . , R • X k spur : set of spurious terms which are selected in the k th run of the search algorithm, but are not present in the actual system, i.e., X k spur = X k run − X • ∅ : the null set On the basis of these definitions, each run of the search algorithm gives one of the following search outcomes: 1. Identification of the Correct Structure (Exact Fitting ) : In this scenario the identified model contains all system terms and does not include any spurious terms, i.e., X k run = X and X spur = ∅

Over Fitting :
The identified model contains all system terms; however spurious terms are also selected, i.e., X k run ⊃ X and X spur = ∅

Under Fitting-1 :
The algorithm fails to identify all system terms; though it does not include any spurious terms, i.e., X k run ⊂ X and X spur = ∅

Under Fitting-2 :
The algorithm fails to identify all system terms; however spurious terms are selected, i.e., X k run ⊃ X and X spur = ∅ It is easier to follow that, qualitatively, the search is successful only when the outcome is Exact Fitting (X alg = X ). A higher number of 'exact-fitting' search outcome indicates a robust search performance.

Test Non-linear Systems
The efficacy of the search algorithms is evaluated through 6 non-linear systems shown in Table 1. These systems have been selected from the existing investigations on the structure detection [9,11,14,27,45,47]. The systems are excited by a white noise sequence having either uniform or Gaussian distribution as shown in Table 1. A total of 1000 input-output data points, (u, y), are generated from each system and 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., 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 (1), i.e., N t = 165.

Results
This study aims to evaluate the performance of Two-dimensional particle swarms (2D-UPSO) on the structure selection problem. For this purpose, the structure selection problem of several non-linear systems described in Section 4.3 is considered. The selection of search parameters and its effects on the search performance of 2D-UPSO are discussed in Section 5.1. For the comparative analysis purposes, four other structure selection algorithms are also considered which include: GA, BPSO, OFR-ERR and OFR-Pruning. The results of this investigation are discussed in Section 5.2.
After establishing the efficacy of 2D-UPSO, it is applied to identify the structure for two wellknown continuous-time nonlinear systems, as discussed in Section 5.4. Finally, in Section 5.5, the 2D-UPSO is applied to real-world application problem of identifying nonlinear wave force dynamics.

2D-UPSO: Selection of Search Parameters
The proposed 2D-UPSO has two tunable search parameters: the unification factor (u f ) and the refresh gap (RG). The unification factor essentially controls the influence of global and local exemplars (i.e., gbest and nbest). As seen in (5), a higher value of u f increases the influence of the global learning set L α 2 and thereby encourages exploitation of the search space. In contrast, a lower u f encourages the exploration of the search space as more weight is assigned to the neighborhood learning set, L α 3 . A detailed discussion on the neighborhood topology of swarm and its influence on the information progression and the search performance can be found in [48].
The refresh gap, RG, determines the interval within which a particle can keep exploring new structures without any improvement over the best structure found hitherto by itself. The refresh gap is introduced as a preventative measure against swarm stagnation, i.e., the situation in which the learning exemplars are trapped in the local minima which can eventually draw the entire swarm into the local minima. To prevent such a scenario, in 2D learning, if a particle cannot determine the better structure for a pre-fixed number of iterations (referred to as Refresh Gap, RG), its velocity is reinitialized, i.e., both cardinality and term selection likelihoods are reinitialized to ensure that the particle can escape from the local minima. The procedure to identify the particle stagnation and velocity re-initialization is outlined in Lines 19-23 and Lines 8-11, Algorithm 1.
To evaluate the effects of u f and RG on the search performance, a grid-search has been carried out using a total of 100 different combinations of {u f , RG} by varying both u f and RG in the following ranges: u f ∈ [0, 1] and RG ∈ [5,50]. For each combination of {u f , RG}, the average (a) S1 + GA (b) S1 + BPSO (c) S1 + 2D-UPSO value of the criterion function, J(· ), over 20 independent runs, is recorded. Note that in this study the structure selection problem is designed as a minimization problem. Hence, a lower value of J(· ) is desirable. For the sake of brevity, only the results obtained with the system S5 and S6 are shown in Fig. 1. The results indicate that the search performance of 2D-UPSO is not very sensitive to the search parameters. Especially, relatively stable search performance is obtained for all test systems in the wide range of u f : [0.3 − 0.7] and RG : [15 − 35]. For the remainder of the study, these parameters are set as u f = 0.4 and RG = 20.

Comparative Evaluation
For each test system in Table 1 determined. For the m th term, x m , the selection frequency is evaluated as follows: selection frequency of the m th term total number of runs It is easier to follow that a higher value of ν is desirable for the system terms. For all test systems, the term selection frequency with all algorithms is determined. However, for the sake of brevity, only the term selection frequency obtained for S1 and S7 is shown in Fig. 2 and 3. Nevertheless, the search behavior for the other systems is not significantly different. Note that at this stage, any spurious terms selected by the algorithms are not removed. Among compared algorithms, it is found that GA is prone to select a higher number of spurious terms with the selection frequency ν >= 0.25, as seen in Fig. 2a-3a. In comparison to GA, the number  of spurious terms selected by BPSO as less, as seen in Fig. 2b-3b. Nevertheless, for most of the systems, both GA and BPSO failed to select some of the system terms. Further, the selection frequency of system terms with GA and BPSO is low,e.g., 0.52 with S2 and 0.67 with S5 (GA); 0.65 with S1 and 0.85 with S5 (BPSO). These results suggest that both GA and BPSO could yield 'under-fitted' structure (as will be discussed later in this section). In contrast, the selection frequency of the system terms by the proposed 2D-UPSO is significantly higher, i.e., ν >= 0.97 for all the systems. Further, 2D-UPSO could yield a significant reduction in the number, and the selection frequency of the spurious terms is obtained as seen in Fig. 2c-3c.
It is worth to emphasize that due to the stochastic nature of the algorithms, the selection of some spurious terms cannot always be ruled out. To investigate if the search performance of GA, BPSO can further be improved, the t-test is applied to identify and remove spurious terms from the term subset/structure obtained in each run of the algorithm, i.e., X k run . Subsequently, the refined structures are used to determine the qualitative search outcome following the definitions in Section 4.2. As discussed earlier, the search for the system structure can be considered successful Nspur 0 0 25 0 0 29 † '' and '' respectively denote inclusion and exclusion of system term; 'Nspur' denotes total number of spurious terms Nspur 1 1 26 0 0 46 † '' and '' respectively denote inclusion and exclusion of system term; 'Nspur' denotes total number of spurious terms only when all system terms have been selected, i.e., when the search outcome is 'Exact-Fitting' (X k run = X , Section 4.2).
The results in Table 2 show the number of times different categories of outcomes, such as 'exact-fitting' or 'over-fitting', are obtained over 40 runs. Note that both GA and BPSO yield under-fitted structures in many runs which indicates the absence of system terms in the selected structure. In contrast, 2D-UPSO could consistently detect true system structure in most of the runs.
In this study, we also compared the performance 2D-UPSO with the classical structure selection approaches such as the Orthogonal Forward Regression (OFR-ERR) [24] and OFR with pruning (OFR-Prune) [45] on all the systems. The limitation of the orthogonal forward regression (OFR) such as the nesting effect, local nature of the search, have extensively been studied in [17,18,[27][28][29]. In contrast, the objective here is to highlight a practical limitation of OFR: the selection of proper threshold, 'σ'. Both, OFR-ERR and OFR-Prune, requires a user-specified σ which determines when the search process will terminate. Hence, the number of terms selected by OFR-ERR is critically dependent on σ. Usually, σ is selected by trial-and-error and in conjunction with the model validity tests [29].
To highlight this issue, the comparative investigation with OFR has been carried out with different values of the threshold. For the sake of brevity, only the results obtained with S5 and S6 are shown in Table 3 and 4. Note that the search outcome is very sensitive to σ and the determination of a proper threshold is not trivial. The improper selection of σ can lead to either exclusion of correct system/significant term or inclusion of many spurious terms.

Performance Under Slow Varying Input
In several existing investigations [9,20,28,29], a worst case scenario for identification has been used as a benchmark to show the effectiveness of their algorithms by considering a system which is excited by a slow-varying input (non-persistent excitation). To evaluate the search performance

Terms
GA BPSO 2D-UPSO Nspur 10 0 0 † '' and '' respectively denote inclusion and exclusion of system term; 'Nspur' denotes total number of spurious terms in such a scenario, the same example has been considered here which is given by, where, 'y(· )' is the observed value of 'w(· )'. The system is excited with a slow-varying input u(· ) which is essentially an AR process with a repeated pole at 0.8. A model set consisting of 165 candidate NARX terms is generated with the following specifications: [n u , n y , n l ] = [4,4,3]. The system is identified by all the compared algorithms following the procedure outlined in Section 4.1. The search outcomes are shown in Table 5. It is observed that both BPSO and 2D-UPSO could identify all the system terms even-though only partial information is extracted from the system due to poor excitation.
It is worth to emphasize that, this example has been considered only for the purpose of comparative analysis and such data should not be used in practice. The system should be persistently excited for identification purposes.

Discrete Models for Nonlinear Continuous-Time Systems
In the previous sections, several discrete polynomial NARX models have been considered to show the effectiveness of the 2D-UPSO. In this section, we apply the proposed approach to fit a discrete time model from the data obtained from continuous time systems. For this purpose, two well-known systems are considered: Duffing's and Van der Pol 's Oscillator. The dynamics of these oscillators are given respectively by,  For both models, the input-output data are generated by exciting the system with a uniformly distributed Gaussian sequence with zero mean and the output is sampled at 500 Hz. The system parameters are set as follows: ω n = 45π, ζ = 0.01 and ε = 3. A total of 1000 data pairs are generated. The discrete NARX model was fitted using 700 data points, and the remaining 300 data points are used for validation.
The initial model is specified with the following NARX parameters: n u = 5, n y = 5, n l = 3, which gives a total of 286 candidate terms. Following the procedure outlined in Algorithm 4, discrete NARX models are identified for both the systems. The identified model of Duffing's and Van der Pole's Oscillator are shown respectively in (12) and (13).

Model Validation using Generalized Frequency Response Functions
The validation of the model is an essential step in the system identification. In all the examples considered in the previous sections the correct structure of the system being investigated was known. Therefore it was straightforward to see whether the correct terms are included or not. However, this is not the case when the data is generated from a continuous time system.
The model is therefore initially validated both by comparing the model predicted output and cross-correlation based model validity tests [49]. In this study, we further establish the accuracy of the models by comparing the Generalized Frequency Response Functions (GFRFs). Note that the discrete model of a continuous time system is not necessarily unique. However, if the discrete model has correctly captured the dynamics of the underlying continuous-time system, then the frequency response functions corresponding to the discrete model should match with those of the continuous time system. This, therefore, has been the motivation to show the efficacy of the proposed approach by comparing GFRFs of the discrete model with the continuous system.
The linear and the third order GFRFs of the continuous time models are shown in Fig. 4a, 4c and 4e (for Duffing's) and in Fig. 5a, 5c and 5e (for Van der Pole's). For the third order GFRF f 3 = f 1 . Note that the GFRF of both the oscillators is characterized by several peaks and ridges as seen in Fig. 4c, 4e, 5c and 5e. The ridges occur both at the linear resonant frequency (±22.5 Hz for both oscillators) and when f 1 + f 2 + f 3 = ±22.5. The peaks occur when all f 1 , f 2 and f 3 are equal to ±22.5 Hz. Further distinguishing characteristic of Van der Pole's oscillator is the presence of 'gorge' along the line given by f 1 + f 2 + f 3 = 0, as seen in Fig. 5c and 5e. It is worth to note that these results are in agreement with the earlier investigation carried out by Billings and Jones [39].
The frequency response of the discrete identified models are shown in Fig. 4b, 4d and   Duffing's) and in Fig. 5b, 5d and 5f (for Van der Pole's). It is clear that both the identified models could accurately capture the frequency response of the continuous models.

Identification of Nonlinear Wave Forces
In this section, a practical problem of non-linear wave force identification is considered. The objective is to investigate if the 2D-UPSO algorithm could identify a practical non-linear system. For this purpose, the velocity and force time histories of a fixed cylinder were obtained from The University of Salford. The input-output data are decimated by a factor of 2 and is shown in Fig. 6. Further information about this data set can be found in [50].
The model is fitted between the in-line force (input, 'u') and the horizontal water particle velocity (output, 'y') considering 900 data points and it is validated with a sequence of over 400 points of input-output data taken arbitrarily from the rest of the available data points. Following similar procedure (Algorithm 3), polynomial NARX model is fitted using GA, BPSO and 2D-UPSO considering the initial NARX model set with [n u , n y , n l , N t ] = [7,7,3,680]. The model is validated by computing the model predicted output over validation data and the cross-correlation plots are shown in Fig. 8. From the results, it is observed that the identified NARX model could successfully encode the dynamics of the wave-force data. Note that although the identified model gave practically acceptable N M SE of 2.12 × 10 −1 , further improvement in the model accuracy may be obtained by fitting a NARMAX model, which includes a noise model, to reduce/eliminate bias in the least square estimates. The proposed 2D-UPSO can be extended to identify NARMAX models which will be reported in the future.

Conclusions
A new Two-dimensional (2D) learning approach of selecting the correct structure of polynomial NARX models has been proposed. One of the key advantage of this algorithm is that the information about the number of terms is explicitly incorporated into the search process. This property significantly reduces the number of spurious terms in the identified model. The effectiveness of the proposed approach has been demonstrated considering both simulated and practical non-linear systems. For discrete-time and practical nonlinear systems, the identified models are validated using the well-known cross-validation and correlation-based model validity tests. For continuoustime nonlinear systems, the identified discrete NARX models are further validated by comparing the generalized frequency response functions (GFRFs) which demonstrate that 2D-UPSO could successfully capture the system dynamics. Let the velocity of i th particle for this problem be given by, This implies that a total of 3 number of terms to be included in the structure. The next step is to determine which 3 terms should be included.