Multi-Objective Evolutionary Framework for Non-linear System Identification: A Comprehensive Investigation

TL;DR

This paper investigates a multi-objective evolutionary framework for selecting structures of nonlinear systems modeled by polynomial NARX models, analyzing the performance of different MOEAs and their robustness across parameters.

Contribution

It introduces a comprehensive multi-objective framework integrating MCDM and evaluates three MOEAs for nonlinear system identification, highlighting their effectiveness and robustness.

Findings

MOEAs can effectively identify nonlinear system structures.

Multiple valid models can be identified for continuous-time systems.

MOEAs are robust over a wide range of parameters.

Abstract

The present study proposes a multi-objective framework for structure selection of nonlinear systems which are represented by polynomial NARX models. This framework integrates the key components of Multi-Criteria Decision Making (MCDM) which include preference handling, Multi-Objective Evolutionary Algorithms (MOEAs) and a posteriori selection. To this end, three well-known MOEAs such as NSGA-II, SPEA-II and MOEA/D are thoroughly investigated to determine if there exists any significant difference in their search performance. The sensitivity of all these MOEAs to various qualitative and quantitative parameters, such as the choice of recombination mechanism, crossover and mutation probabilities, is also studied. These issues are critically analyzed considering seven discrete-time and a continuous-time benchmark nonlinear system as well as a practical case study of non-linear wave-force modeling. The results of this investigation demonstrate that MOEAs can be tailored to determine the correct structure of nonlinear systems. Further, it has been established through frequency domain analysis that it is possible to identify multiple valid discrete-time models for continuous-time systems. A rigorous statistical analysis of MOEAs via performance sweet spots in the parameter space convincingly demonstrates that these algorithms are robust over a wide range of control parameters.

Tables

Table 1. Table 8: Performance Metric ( P ​ M 𝑃 𝑀 PM ) for System 𝒮 1 subscript 𝒮 1 \mathcal{S}_{1}

Parameter Combination $\{{𝒑}_{𝒄}^{𝒋},{𝒑}_{𝒎}^{𝒋}\}$	NSGA-II		SPEA-II		MOEA/D
	Single Point	Uniform	Single Point	Uniform	Single Point	Uniform
$\{0.1,0.001\}$	0.9696	0.9764	0.9649	0.9419	0.5240	0.5573
$\{0.1,0.0035\}$	0.9836	0.9867	0.9774	0.9883	0.9830	0.9820
⋮	⋮	⋮	⋮	⋮	⋮	⋮
$\{0.8,0.001\}$	0.9812	0.9890	0.9542	0.9704	0.8668	0.9203
⋮	⋮	⋮	⋮	⋮	⋮	⋮

Table 2. Table 9: Statistical Comparison of Crossover Mechanisms

System	Null Hypothesis (${𝑯}_{\mathrm{𝟎}}$)
	NSGA-II	SPEA-II	MOEA/D
${𝒮}_1$	✗	✓	✗
${𝒮}_2$	✓	✓	✗
${𝒮}_3$	✓	✗	✗
${𝒮}_4$	✓	✓	✗
${𝒮}_5$	✗	✗	✗
${𝒮}_6$	✗	✓	✗
${𝒮}_7$	✗	✗	✗

Equations

$$y(k)$$

$$n$$

$$y(k)$$

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

$${\vec{J}}_1=\{J_1({\mathcal{X}}_1),J_2({\mathcal{X}}_1),\dots J_m({\mathcal{X}}_1)\}{\vec{J}}_2=\{J_1({\mathcal{X}}_2),J_2({\mathcal{X}}_2),\dots J_m({\mathcal{X}}_2)\}$$

$$\displaystyle if{f}\ \$$

$${\mathcal{X}}_{model}=\left\{\begin{array}{cccc}{x}_1& x_2& \dots & x_n\end{array}\right\}$$

$${\mathcal{X}}^{\star }$$

$$\vec{J}({\mathcal{X}}_i)$$

$$\text{where, }{\xi }_i$$

$$\displaystyle\Gamma^{\star}=\Big{\{}\mathcal{X}\subset\mathcal{X}_{model}\ |\ \nexists{\mathcal{X}_{j}}\subset\mathcal{X}_{model}\mathrel{\mathop{\mathchar 58\relax}}\mathcal{X}_{j}\preceq\mathcal{X}\Big{\}}$$

$${\beta }_i$$

$${\mathcal{X}}_{model}$$

$$=\left\{\begin{array}{ccccc}y(k-1)& u(k-2)& y(k-3)& y(k-2)u(k-2)& u(k-3{)}^3\end{array}\right\}$$

$${\beta }_i$$

$${\mathcal{X}}_i$$

$$\vec{J}({\mathcal{X}}_i)$$

$$\text{with, }\mathcal{P}$$

$${\mathcal{S}}_1\text{mathchar}58y(k)=$$

$$-0.7y(k-2)u(k-2{)}^2-0.5y(k-2)+e(k)$$

$$a_{i,j}$$

$$\text{where, }{\delta }_O$$

$$\vec{w}$$

$$\left[\begin{array}{cc}{a}_{\xi ,\xi }& a_{\xi ,\mathcal{E}}\\ a_{\mathcal{E},\xi }& a_{\mathcal{E},\mathcal{E}}\end{array}\right]=\left[\begin{array}{cc}1& 5\\ \frac{1}{5}& 1\end{array}\right]$$

$$w_{\xi }$$

$${\vec{J}}^{\star }$$

$$\mathcal{D}({\mathcal{X}}_i)$$

$$T_p({\mathcal{X}}_i,\Gamma )$$

$$\text{where,}{t}_p({\mathcal{X}}_i,{\mathcal{X}}_j)$$

$$\displaystyle\mathcal{R}(\mathcal{X}_{i})=\Big{(}\prod\limits_{p=1}^{m}T_{p}(\mathcal{X}_{i},\Gamma)^{w_{p}}\Big{)}^{\frac{1}{m}}$$

$${\Gamma }_A$$

$${\Gamma }_B$$

$$\mathcal{C}(A,B)$$

$$\vec{r}$$

Full text

††journal: Neurocomputing

\DeclareFloatFont

henry\floatsetup[figure]font=henry

1 Introduction

The system identification aims to encapsulate dynamical properties of the system under scrutiny in a mathematical descriptor/model. Such identified models are not only critical to understand the system dynamics but also to the subsequent prediction and control applications. Given that the dynamic input-output behavior of most of the systems is predominantly nonlinear in nature, the identification of nonlinear systems has attracted significant attention over the past few decades [1, 2, 3, 4, 5, 6, 7, 8]. The first step of the identification process is the selection of appropriate system representation, e.g., Volterra series, Wiener, Hammerstein, Neural Networks, polynomial Nonlinear Auto-Regressive Moving Average with eXogenous inputs (NARMAX) and others. Among these representations, the NARX models [9, 10, 4, 11] are arguably the most popular due to various advantages, e.g., convenient linear-in-parameter form and the availability of frequency domain analysis tools [12, 4]. The identification of non-linear systems in the form of a polynomial NARX models is, therefore, considered in this investigation.

The system identification is essentially a data-driven modelling process which begins by defining a search space of candidate terms/basis functions. The subsequent model building procedure involves mainly two steps: 1) Structure Selection: selection of significant terms/basis functions and 2) Parameter Estimation: estimation of associated parameters. Due to linear-in-parameter form of the NARX model, the parameters can be estimated using the least-squares based approaches, provided the terms to be included in the model are accurately determined. However, the selection of correct terms which captures the system dynamics is one of the major challenges of the system identification [6, 7], and, therefore, is the main focus of this investigation.

The exhaustive search to solve the structure selection problem is intractable even for moderate number of terms ‘$n$’, as this would require an evaluation of $2^{n}$ term subsets/structures. The structure selection problem is often ‘NP-Hard’ [13] and therefore requires an efficient search strategy. Since the search for system structure is primarily dependent on the limited number of observed data, it involves bias-variance trade-off [2, 3, 7], i.e., too sparse a structure may yield a high bias error (under-fitting) whereas a very complex structure may yield a high variance error (over-fitting). Further, earlier investigations [14, 15] show that the over-fitted structures tend to induce ‘ghost’ dynamics which are not present in the actual system. Pursuing these arguments, it is clear that determination of the size of the identified structure (i.e., number of terms or ‘cardinality’) is a critical aspect of the structure selection process. However, this issue has received relatively less attention in the existing investigations.

Among the existing structure selection approaches, the Orthogonal Forward Regression (OFR) proposed by Billings and Korenberg [16] has extensively been studied. It is a sequential model building approach in which the most significant term is added in each step, which is determined on the basis of a statistical criterion. Over the years, this notion of the original OFR approach has been refined to further improve the search performance. These include OFR approaches with either new term evaluation criteria [17, 18, 19, 20] or improved search strategies [21, 22, 23, 24, 25]. A detailed treatment on this subject can be found in the recent investigation by the authors [25]. It is worth noting that, while OFR and its extensions are quite effective in detecting the significant terms, these approaches typically require an auxiliary procedure to estimate the number of terms (cardinality). Usually, an information theoretic criterion is used for this purpose, e.g., see [26]. However, it is known that the core assumptions of the information criteria may not always hold. In such scenarios, the information criteria usually tend to over-fit as suggested by the earlier investigation of Nakamura et al. [27]. Hence, the estimation of cardinality still requires further investigations in OFR based approaches.

The other interesting approach to structure selection is based on stochastic sampling of the search space, e.g., Evolutionary Algorithms [28, 29, 30, 31, 32, 33, 34] and Bayesian inference [35, 36]. In the earlier investigations, Genetic Algorithm (GA) [28, 29] and Genetic Programming (GP) [30, 31] have been proposed as an alternative to OFR. In these investigations, the structure selection has been posed as a single objective problem, where the search objective is chosen to be an aggregated function of predictive performance and number of terms to balance the bias-variance dilemma. Recently, a randomized structure selection approach is proposed in [32, 33, 34]. In this approach, each candidate term is assigned with selection probability which is updated following a positive-feedback mechanism similar to ant colony optimization [37]. In [36], an alternative approach based on Bayesian framework is proposed wherein the posterior probabilities of both the terms and parameters are determined through Markov Chain Monte Carlo sampling.

Although the existing approaches have been shown to be effective, it is clear that their primary search focus is to determine the significance of the terms. The pursuit for parsimonious structure with a better predictive performance is essentially contradictory in nature. Therefore, the most challenging task of the structure selection is, arguably, to obtain a good balance of bias-variance dilemma. It is known that a unique solution to such problems may not exist. For instance, a discrete NARX model of a continuous-time system may not be unique [4]. Pursuing these arguments, it is clear that the structure selection is, in essence, a multi-objective problem. Note that over the past two decades, several Multi-Objective Evolutionary Algorithms (MOEAs) have specifically been developed to address such problems, e.g., see [38, 39, 40]. Hence, in this scenario, the following question arises: can the existing MOEAs be used to reach a better trade-off in bias-variance dilemma?

The present study aims to answer this pertinent question. In particular, this study proposes a Multi-Objective Structure Selection (MOSS) framework which integrates all the key components involved in the Multi-Criterion Decision Making (MCDM) process to address both aspects of structure selection, i.e., term significance and cardinality estimation.

It is worth noting that formulation of the structure selection problem as a multi-objective problem is not new, but it is under-explored. The research in this direction was pioneered in the early nineties by Fonseca et al. [28, 41, 42, 43]. In [28, 41], the structure cardinality is pre-specified and the structure is identified to balance multiple prediction criteria. Subsequently, a multi-objective GP is developed for the structure selection on the basis of classical Multi-objective GA (MOGA) [44] which is extensively investigated in [42, 43, 45, 46]. Since then, several seminal MOEAs have been developed with significant improvement in the search performance (e.g., see [38, 39, 40]). However, in recent years, the application of MOEAs in structure selection has received relatively less attention except for a few preliminary investigations in [47, 48, 49]. While the existing investigations have demonstrated the capabilities of MOEAs to address the structure selection problem and provided proof-of-concept, several challenging issues are yet to be investigated. The focus of this study is, therefore, to address some of these issues,

• •

Investigate the performance of well-known MOEAs such as NSGA-II, SPEA-II and MOEA/D, and to determine whether there exists any significant difference through rigorous statistical comparison.

• •

Determine the sensitivity of MOEAs to various qualitative and quantitative control parameters and determine the robustness through performance sweet spots [50, 51] in the parameter space.

• •

Evaluate the a posteriori MCDM techniques through Generalized Frequency Response Functions (GFRFs) in frequency domain validation.

For this purpose, the performance of MOEAs is critically analyzed from the identification perspective considering seven discrete-time and a continuous-time benchmark system as well as a practical identification problem in structural hydro-dynamics.

It is worth noting that this study is a part of our broader investigation on alternative approaches for structure selection. In the first part of this investigation [25], the orthogonal floating search algorithms have been proposed to alleviate the nesting effects of the classical OFR. Subsequently, a new two-dimensional structure selection approach was developed in [52, 53], which explicitly integrates the information about subset cardinality (number of terms) into a probabilistic learning framework to balance the bias-variance dilemma. These approaches have been shown to be very effective on both simulated and practical nonlinear systems. These algorithms, however, have been designed to identify a unique structure which is optimum for certain criterion. The multi-objective framework (MOSS) proposed in this study can alleviate this issue. MOSS can identify multiple valid models for the given systems and thereby yield a better trade-off for the bias-variance dilemma. This, not only aids the decision making process but also instills confidence in the identified models.

The rest of the article is organized as follows: the polynomial NARX model and the concept of Pareto dominance are discussed briefly in Section 2. This is followed by Multi-objective Structure Selection (MOSS) framework in Section 3, which describes, in detail, all the components involved in the multi-criteria decision making (MCDM) process. The framework of this investigation is introduced in Section 4 which provides the detail about the benchmark nonlinear test systems, search environment and the performance metrics being used in the comparative evaluation. The results of this investigation are discussed in three stages. First, the MOEAs are critically analyzed from the identification perspective on the benchmark systems in Section 5. Next, the search performance of MOEAs is compared through rigorous statistical analysis in Section 6. Finally, the robustness of MOEAs to the search parameters is analyzed at length in Section 7, followed by the conclusions in Section 8.

2 Preliminaries

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 Pareto Dominance

For multi-objective problems, the ideal solution which would yield optimum values across all the objectives may not exist; rather there may exist a set of non-dominated or Pareto Optimal solutions [54]. Such non-dominated solutions may not necessarily be optimum for each objective, however, they are better than the other solution when all the objectives are simultaneously considered.

To understand the concept of the Pareto dominance, consider a problem with ‘$n$’ decision variables and having a total of ‘$m$’ number of objectives. Let $\mathcal{X}_{1}$ and $\mathcal{X}_{2}$ denote the candidate solutions of this problem. Further, denote the corresponding objective function vectors by $\vec{J}_{1}$ and $\vec{J_{2}}$, as follows:

[TABLE]

The Pareto dominance of these solutions can be determined on the basis of the objective values as follows: $\mathcal{X}_{1}$ dominates $\mathcal{X}_{2}$ (denoted as $\mathcal{X}_{1}\preceq\mathcal{X}_{2}$),

[TABLE]

In this study, the set of such non-dominated solutions is denoted by ‘$\Gamma$’ and the corresponding set of objective function vectors by ‘$\Lambda$’.

3 Multi-objective Structure Selection Framework

In this study, the structure selection problem is posed as a Multi-Criteria Decision Making (MCDM) problem to balance the bias-variance dilemma. The objective here is to minimize both the predictive error as well as the number of terms (cardinality). Without loss of generality, this study takes an a posteriori approach to the multi-criteria decision making, where a set of non-dominated structures are first identified by an MOEA. Next, an a posteriori decision method is used to select the final structure from the pool of identified structures. Such approach decouples the search and decision making process. Since the search process is unbiased to any of the objective, it is possible to accommodate distinct preferences a posteriori without the need to reiterate the search process.

To this end, a Multi-objective Structure Selection (MOSS) framework is proposed which consists mainly of three procedures: 1) Preference Articulation 2) Search for Non-dominated Structures and 3) a posteriori MCDM.

The search objectives of the proposed MOSS are discussed in Section 3.1. To solve this problem, three seminal Multi-objective Evolutionary Algorithms (MOEAs) are being considered: NSGA-II [38], SPEA-II [39] and MOEA/D [40]. These are briefly discussed in Section 3.3. In all the MOEAs, the binary string representation is used to encode system structures, as discussed in Section 3.2. Further, in system identification, the Decision Maker (DM) is often interested in a particular region of the Pareto front. This information is incorporated to focus the search in the region of interest by introducing a priori goal point, which is discussed in Section 3.4. The search process yields a set of non-dominated structures, which is referred to as Approximate Pareto Set (APS). For the final selection, the identified non-dominated structures in APS are ranked using either Minimum Manhattan Distance (MMD) or Multi-criteria Tournament Decision (MTD) approach, which is discussed at length in Section 3.5.

All the steps involved in the proposed MOSS are outlined in Fig. 1. In the following, these steps are discussed in detail from the identification perspective.

3.1 Multi-objective Structure Selection Problem

To understand the structure selection problem, consider the identification of a system from a total ‘$\mathcal{N}$’ number of input-output data-points, i.e., $\{u_{i},y_{i}\},i=1,2,\dots,\mathcal{N}$. Denote the given super-set of NARX terms by ‘$\mathcal{X}_{model}$’, as follows:

[TABLE]

where, $x_{i}\in\mathcal{X}_{model}$ could represent any linear or non-linear term of the NARX model in (2); ‘$n$’ denotes the total number of NARX terms which is expressed in  (1). It is dependent on the maximum lag of inputs and outputs as well as on the degree of non-linearity.

For the given identification dataset and a model-set $\mathcal{X}_{model}$, the goal of the structure selection is to identify the subset of system/significant terms, $\mathcal{X}^{\star}\subset\mathcal{X}_{model}$. This can be formulated as an optimization problem as:

[TABLE]

where, ‘$J(\cdot)$’ denotes a suitable criterion function. Most of the existing structure selection approaches often use the predictive performance of the model as the criterion function [28, 29, 30, 31, 32, 33, 34]. These can, therefore, be categorized as the single objective approach to structure selection. In contrast, in this study, the structure selection is formulated as a multi-objective optimization problem with the following search objectives:

[TABLE]

[TABLE]

In (7), ‘$\mathcal{X}_{i}$’ denotes the $i^{th}$ structure under consideration; ‘$\xi_{i}$’ denotes the cardinality (number of terms); ‘$\mathcal{E}_{i}$’ denotes the Normalized Mean Squared Error (NMSE) obtained over the validation data; ‘$\hat{y}$’ denotes the model predicted output; and ‘$\mathcal{N}_{v}$’ denotes the length of the validation data.

It is easier to follow that the search objectives, $\{\xi,\ \mathcal{E}\}$, have been selected to identify parsimonious models which would give better predictive performance. For such multi-objective formulation, there may exist several non-dominated structures instead of single optimum structure ($\mathcal{X}^{\star}$) (as discussed in Section 2.2). Hence, the goal of the multi-objective structure selection is to identify the set of non-dominated structures, $\Gamma^{\star}$, which are Pareto optimal for the criteria function, $\vec{J}(\cdot)$ in (6):

[TABLE]

It is worth noting that the formulation of the structure selection as a multi-objective optimization problem is not new. The pioneer investigation in this direction has been carried out by Fonseca et al. in [28, 41, 43, 45]. In most of these investigations, the structure selection problem is essentially formulated as a many-objective problem with $5$ or more search objectives, i.e., $m\geq 5$. For instance, in [45], to reduce the model complexity, maximum lags and degree of nonlinearity are also included as the search objectives. Further, linear and higher order cross-correlation functions between input and error residuals are also framed as the search objectives to embed the model validity tests [55] into the search process. While the inclusion of these additional objectives can narrow the feasible search space, in this study, we have limited to number of search objectives to two (i.e., $m=2$) for the following reasons:

• •

The identification of the Pareto optimal structures is essentially the penultimate step. In the last step, the decision maker (DM) has to select the final structure from the pool of non-dominated structures which is often referred to as a posteriori multi-criteria decision making (MCDM). By limiting the search objectives to $\{\xi,\ \mathcal{E}\}$, a posteriori MCDM can be reduced to a model order selection problem.

• •

The dominance relations tend to lose the comparative power for a higher number of objectives. This could deteriorate the search performance of the dominance based MOEAs such as NSGA-II and SPEA-II (which are being used in this study) [54, 51].

• •

The complexity of preference articulation is directly dependent on the number of objectives. For instance, the specification of multiplicative preference relations to select the final structure becomes increasingly difficult as the number of objectives increase.

3.2 Solution Representation

Let the set containing all model terms be denoted as $\mathcal{X}_{model}=\begin{Bmatrix}x_{1}&x_{2}&\dots&x_{n}\end{Bmatrix}$. For this problem, each member of an MOEA ($\beta$) encodes a candidate structure or set of terms in an $n$ - dimensional binary vector. For example, the $i^{th}$ member, $\beta_{i}$, is represented as:

[TABLE]

The $m^{th}$ term ($x_{m}$) from $\mathcal{X}_{model}$ is included into the candidate structure provided the corresponding bit in the member, ‘$\beta_{i,m}$’ is set to ‘$1$’.

For sake of clarity, consider a simple NARX model with a total of $5$ terms ($n=5$) as follows:

[TABLE]

For this problem, assume that the $i^{th}$ member is given by,

[TABLE]

This implies that only the first, third and fourth terms from the set $\mathcal{X}_{model}$ are included into the structure/term subset. Thus, the structure ‘$\mathcal{X}_{i}$’ encoded by the particle $\beta_{i}$ is given by,

[TABLE]

3.3 Multi-objective Evolutionary Algorithms

The Multi-Objective Evolutionary Algorithms (MOEAs) are designed to address the multi-criteria decision making (MCDM) process a posteriori. It is, therefore, vital that the search process identifies well distributed solutions which are closer to the true Pareto front. This provides a better insight to the Decision Maker (DM). These two main goals of the search process are, however, contradictory in nature. The existing MOEAs can be categorized based on their approach to balance the convergence with the solution diversity of the approximated Pareto front. In this study, the following three seminal MOEAs are selected to represent two distinct categories of MOEAs: Pareto dominance based (NSGA-II [38] and SPEA-II [39]) and decomposition based (MOEA/D [40]).

In the following subsections, first the solution representation used to encode term subset/structure is discussed. Next, the compared MOEAs are briefly summarized.

3.3.1 Pareto Dominance based MOEA

Most of the existing Pareto Dominance based MOEAs are based on Goldberg’s [56] earlier proposal of the direct integration of dominance relations into the search process. This integration can loosely be thought of as a pseudo fitness approach, where a scalar value derived from the dominance relations is assigned to the members to accommodate multiple search objectives. However, since the Pareto optimality defines only a partial order [54, 57], several members of the population may become incomparable due to similar value of the pseudo fitness, e.g., all the non-dominated members are incomparable. Such a scenario may lead to genetic drift phenomenon [56], i.e., even though all the non-dominated members are equivalent, the population tend to converge only to the part of the Pareto front. An additional distinguishing metric is, therefore, required. To this end, different measures are introduced which essentially encapsulates the diversity information, e.g., fitness sharing, crowding distance, density estimate. The idea here is to assign more importance to the members in an under-represented part of the Pareto front and thereby promote the solution diversity. Furthermore, it has been shown that the introduction of elitism can lead to improved convergence properties [58].

In this study, the following two seminal dominance based MOEAs are considered: Non-dominated Sorting Genetic Algorithm-II (NSGA-II) [38] and Strength Pareto Evolutionary Algorithm-2 (SPEA-II) [39]. While both NSGA-II and SPEA-II are based on the Goldberg’s proposal, they differ in terms of how these search goals are achieved.

NSGA-II can be thought of as a direct extension of the Goldberg’s proposal. In NSGA-II, the dominance relations are used to rank and segregate the entire population into successive fronts using a well-known non-dominated sorting operation. Each front contains the members with a similar dominance rank. To promote the solution diversity, a crowding-distance metric is proposed which essentially determines the average distance to neighbors in the objective space. The mating pool is determined using the Crowded Tournament Selection (CTS) which is similar to binary tournaments except with a modified selection operator to balance convergence with diversity, as follows. The member with a lower rank is preferred, however, when the ranks are incomparable, the tie is broken by selecting a member with smaller crowding distance. The offspring population is generated by usual recombination and mutation operators. The elitism is maintained by combining population of ‘$ps$’ parents with the population of ‘$ps$’ offspring to form the population union with $(2\times ps)$ members. The elitism is ensured by allowing only ‘$ps$’ members of the union population to survive where the non-dominated members are given first priority followed by other members with the lower rank and smaller crowding distance.

In contrast, SPEA-II maintains elitism by keeping an external population (referred to as archive) of non-dominated solutions which is continually updated throughout the search process. Although, the archive size may differ from the population size, in this study both are fixed to the same value, and is denoted by ‘$ps$’, for the sake of simplicity. Further, in SPEA-II, dominance relations are incorporated by introducing the concept of ‘strength value’. For each member in the population and in the archive, the strength value is determined from the total number of dominated members. Further, the Euclidean distance in objective space to the $k^{th}$ nearest neighbor is used as the density estimate, where ‘$k$’ is the function of the population size and the archive size. To balance convergence with diversity, the pseudo fitness in SPEA-II is composed of two components: ‘raw’ fitness derived from strength value and the density estimate. The mating pool is generated by the binary tournament selection among the members of the archive. The offspring population is generated by usual recombination and mutation operators. Consequently, the external archive is updated to maintain elitism. To this end, in each iteration, the archive and the offspring population are combined to form the union population of size $(2\times ps)$. The elitism is maintained by including all the non-dominated members from the population union into the archive. Note that in SPEA-II, the archive size is fixed i.e., the archive is always filled. Accordingly, two different scenarios may exist based on the number of non-dominated solutions. If the size of the archive, $ps$, is higher than the number of non-dominated members, the remaining positions in the archive are filled by selecting the members with the lowest pseudo-fitness from the population union. On the other hand, if the non-dominated members are more than ‘$ps$’, then the additional members from the crowded region of the Pareto front are removed using the density estimate metric.

Since NSGA-II and SPEA-II have extensively been studied from both theoretical and empirical perspective, further details and pseudo codes for these algorithms are omitted here, and can be found elsewhere [38, 54, 39]. Nevertheless, in the following the necessary implementation details are discussed. Since one of the objectives of this study is to compare the search performance of MOEAs, the reproduction is performed using similar recombination and mutation operators in all the compared MOEAs. Hence, the first step is to select the appropriate reproduction operators for the binary string representation which is used to encode the candidate term subsets/structures (see Section 3.2). To this end, two recombination operators are considered in this study: single-point and uniform crossover. A detailed investigation to study the effects of recombination operator will be discussed in Section 7.3. Based on this investigation, the uniform crossover has been selected in this study, unless specified otherwise. Further, the mutation is induced through flip-bit mutation.

The common reproduction operators used in this study for NSGA-II and SPEA-II are shown in Algorithm 1, where ‘$\beta_{i}$’ and ‘$\hat{\beta}_{i}$’ respectively denote the $i^{th}$ parent and the corresponding offspring. Note that the selection procedure in NSGA-II and SPEA-II are different. NSGA-II uses Crowded Tournament Selection (CTS) which is based on ranks and the crowding distance as discussed earlier in this section. In contrast, SPEA-II uses Binary Tournament Selection (BTS) on the basis of the pseudo fitness function. Further issues regarding the parameter settings such as crossover and mutation probabilities will be discussed at length in Section 7.2.

3.3.2 Decomposition based MOEA

MOEA/D takes a different approach to balance convergence with diversity. In MOEA/D, each member of the search population is tasked with solving a unique scalar sub-problem of a decomposed multi-objective problem. To this end, a multi-objective problem is decomposed into several scalar sub-problems by a set of uniformly spaced weight vectors, ‘$\lambda$’. For instance, with a population of ‘$ps$’ members, an $m$-dimensional problem is decomposed into a total of ‘$ps$’ sub-problems. Each sub-problem corresponds to a unique weight vector, e.g., the $i^{th}$ weight vector, $\lambda^{i}$, corresponds to $i^{th}$ sub-problem, where, $\lambda^{i}=\{\lambda_{1}^{i},\ \lambda_{2}^{i},\ \dots,\ \lambda_{m}^{i}\}$ and $\sum\limits_{p=1}^{m}\lambda_{p}^{i}=1$. The candidate solution for each sub-problem is evaluated using the corresponding weight vector and an aggregation function such as weighted sum function or Tchebycheff function [40].

Further, prior to the search process, for each sub-problem, a neighborhood of ‘$T$’ sub-problems is determined such that the corresponding distance in the weight vectors is minimized. The core idea of MOEA/D is to exploit the neighborhood relations for the benefit of the search process, i.e., if the $j^{th}$ sub-problem is in the neighborhood of the $i^{th}$ sub-problem, then a good candidate solution for the $i^{th}$ sub-problem may be promising for the $j^{th}$ sub-problem as well. In each iteration, a new candidate solution is determined by reproduction operations on two randomly selected members of the same neighborhood. The new solution is allowed to replace a total of ‘$nr$’ ($nr<T$) worse solutions in the neighborhood. Note that the population in MOEA/D contains the best solution found hitherto for each sub-problem. In addition, an external archive of non-dominated solutions is also maintained and updated in each iteration.

For the sake of fair comparison, the recombination and the mutation operator being used with MOEA/D are similar to NSGA-II and SPEA-II, and are shown in Algorithm 1. However, the selection procedure in MOEA/D is different, where two parents are selected at random from the neighborhood of the parent under consideration.

3.4 Preference Articulation

The preferences of the Decision Maker (DM), if incorporated correctly, can simplify the MCDM process. To this end, the preferences can be formulated a priori, a posteriori or progressively with respect to the search process [44, 59, 60, 61, 62, 63, 64]. This study takes a two-pronged approach to the preference articulation wherein the a priori preferences are formulated as goals in the objective space [44, 59] and the a posteriori preferences are formulated as weights [65]. This distinction arises from the role being played by the preferences with respect to the search process. The key role of the a priori preferences is to focus on the subset of actual Pareto front to better suit the needs of the system identification (to be discussed in Section 3.4.1). In contrast, the a posteriori preferences are being used to select the final structure from the set of identified non-dominated structures (to be discussed in Section 3.4.2 and 3.5.2).

It is worth emphasizing that the preference articulation is not critical to the search performance of the MOEAs. Instead, they are crucial to simplify the a posteriori MCDM process by focusing on a particular region of the actual Pareto front which is of interest to the DM.

3.4.1 A priori Preferences

Most of the existing MOEAs are designed to support MCDM in an a posteriori manner [38, 54, 39, 40], i.e., the goal is to balance diversity in the non-dominated solutions with the convergence to the optimal Pareto front. For instance, in NSGA-II [38], the crowding distance operator has specifically been designed to encourage diversity in the identified non-dominated solutions. This is primarily based on the assumption that the DM does not have any a priori preferences.

However, this may not be the case in system identification applications, wherein it is possible to formulate partial a priori preferences. For instance, earlier investigation suggests that most of the non-linear dynamics can sufficiently be described by $20$ or fewer NARX terms [11], i.e., $\xi\leq 20$. Further, it is easier to follow that a similar upper limit can be determined for the other objective, i.e., NMSE ($\mathcal{E}$). In practice, the structures with NMSE higher than $30\%$ are usually not recommended, as this indicates a high bias error [3]. This, in essence, specifies the region of interest for most of the identification applications.

If such information can be embedded in the search process, then MOEAs can focus only on the portion of the Pareto front which is of interest to the Decision Maker (DM). While this exercise has several benefits (e.g., speed, focus) [60], in this study particular interest is on ‘focus’. By focusing the search only to regions of interest to the DM, irrelevant structures in the approximate Pareto front can be reduced which can simplify the a posteriori MCDM process. This has been the main motivation for including preference articulation in this study.

The integration of DM preferences in MOEA has been investigated in several existing investigation, e.g., see [44, 59, 60, 61, 62, 63, 64]. The preferences can be formulated as weights, priorities or goals [59, 61, 60]. The weights and priority essentially assign the importance to each objective as per preference of the DM. In goal formulation, the preferences are specified in the objective space as an utopian or minimum levels to be attained by the objectives. For the system identification, articulation of a priori preferences as goals is more appropriate for the following reasons:

• •

Specification of goal in objective space is relatively easier than determining the priority or weight for each objective a priori.

• •

The approximate Pareto set obtained with the goal formulation is essentially subset of the original Pareto set, and therefore not biased towards any objective.

Hence, without the loss of generality, it is possible to formulate a priori preferences for the structure selection problem as the following goals: $\xi\leq 20$ and $\mathcal{E}\leq 30$, i.e., the ‘goal’ point for the structure selection problem can be given by: $\{\xi_{lim},\mathcal{E}_{lim}\}=\{20,30\}$, where ‘$\xi_{lim}$’ and ‘$\mathcal{E}_{lim}$’ denote respectively the maximum allowable limit for the cardinality and the prediction error.

There exist several approaches to focus the search efforts in the specific region which is determined by the specified goal. For example, Fonseca and Fleming [44, 59] proposed new dominance relations wherein a higher priority is assigned to the objective in which the goal is not attained. While this approach can easily be integrated with the dominance based MOEAs such as NSGA-II and SPEA-II, it may be challenging for non-dominance based MOEAs such as MOEA/D. Hence, in this study the other approach to goal attainment is considered which casts goals as constraints in the objective space as follows:

[TABLE]

where, the bracket-operator, ‘$\langle\cdot\rangle$’ in (14) returns the absolute value if the operand is negative else it returns zero. The penalty function, $\mathcal{P}$, has been designed to discourage exploration beyond the specified goal point, $\{\xi_{lim},\mathcal{E}_{lim}\}$.

The effects of a priori preferences is illustrated through the following numerical example:

[TABLE]

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

To investigate the effects of goal formulation, NSGA-II is applied to obtain the approximate Pareto fronts in two scenarios: 1) Without any a priori preferences, i.e., using the criteria function in (6) and 2) With a priori preferences, i.e., using the criteria function in (13). For the second scenario, the goal point is specified as $\{\xi_{lim},\mathcal{E}_{lim}\}=\{20,30\}$. The approximate Pareto front obtained over $40$ independent runs of NSGA-II for both the scenarios are shown in Fig 2. It is clear that with goal formulation, the search is focused in the region of interest. Due to the penalty function, $\mathcal{P}$, the structure outside the region of interest ‘appear’ less rewarding. Hence, the search is focused in the boundaries defined by the specified goal point. Further, note that the front identified with the goal formulation is the subset of the approximate Pareto front which is obtained without preference articulation, i.e., no bias is introduced towards any of the search objective.

This procedure of accommodating the goal formulation into the objective space is outlined Algorithm 2.

It is worth noting that, without the loss of generality, the goal point $\{\xi_{lim},\mathcal{E}_{lim}\}$ is fixed to $\{20,30\}$ on the basis of prevailing common knowledge in system identification. While this can be considered as a rule of thumb, some exception may exist. A simple approach to identify such exceptions and the remedial procedure are discussed in Section 5.2.

3.4.2 A Posteriori Preferences

Once the search process is complete and the approximate Pareto set is obtained, the next step is to select the system structure from the pool of non-dominated structures, as shown in Fig 1. To this end, the DM may be impartial towards the objectives or may have particular preferences. The selection of an appropriate a posteriori MCDM approach is therefore dependent on the prevailing scenario, and will be discussed in detail in Section 3.5. Further, the human preferences are often partial and/or abstract [62, 60, 65]. Hence, for the scenario where the DM is biased towards a particular objective, such partial preferences have to be quantified before they can be used with a suitable MCDM approach. While there exist several approaches for this purpose (e.g., see [61, 60, 65]), in this study, a posteriori preferences are formulated as weights using the preference ordering approach [65, 66] due to its simplicity. In the following, this approach is discussed briefly in the context of system identification.

In the preference ordering approach [65, 66], as the name suggests, each objective is ranked based on the order of its importance by the DM. For instance, the DM may prefer parsimonious structures over predictive performance, i.e., ‘$\xi$’ is preferred over ‘$\mathcal{E}$’. Consequently, the objectives are ranked as follows: $[O_{\xi},O_{\mathcal{E}}]=[1,2]$, where ‘$O$’ denotes the objective rank. It is easier to follow that, for a total of $m$ number of objectives, an $m$-dimensional ranking vector is required.

Next, the intensity of the preference, denoted by $\mathcal{I}$, is fixed on a scale from ‘$1$’ to ‘$9$’; where $\mathcal{I}=1$ denotes equal importance to the objectives and $\mathcal{I}=9$ denotes extreme prejudice. Subsequently, the multiplicative preference relations (denoted by ‘$a$’) are obtained using the objective rankings ($O$) and the preference intensity ($\mathcal{I}$) as follows:

[TABLE]

and ‘$a_{i,j}$’ denotes the multiplicative preference over objectives $i$ and $j$. This implies that the $i^{th}$ objective is $a_{i,j}$ times more important than the $j^{th}$ objective.

Finally, the preference weights (denoted by $w$) are determined from the aggregated multiplicative preferences:

[TABLE]

For illustration purposes, consider a particular identification problem. Let the objective ranks and the preference intensity specified by the DM be given by: $[O_{\xi},O_{\mathcal{E}}]=[1,2]$ and $\mathcal{I}=5$. Using (16), the multiplicative preference relations can be determined as:

[TABLE]

Consequently, the preference weights determined using (17) are as follows:

[TABLE]

The complexity of this procedure increases with the number of objectives. It is worth emphasizing that, by limiting number of objectives to two (i.e., $m=2$), the determination of preference weights is significantly simplified. This, in part, is the motivation behind limiting the of objectives to cardinality ($\xi$) and prediction error ($\mathcal{E}$).

3.5 A Posteriori Multi-Criteria Decision Making

Given a set of non-dominated structures, the objective of the a posteriori Multi-Criteria Decision Making (MCDM) methods is to rank these structures to facilitate the final selection by a Decision Maker (DM). Note that the DM may or may not have particular preferences. To accommodate these two distinct scenarios, the following two MCDM methods have been selected in this study: 1) Minimum Manhattan Distance (MMD) [67] and 2) Multi-criteria Tournament Decision (MTD) [66].

Note that both MMD and MTD are independent of the search algorithm. Hence, they can easily be combined with any MOEA for a posteriori decision making. In the following subsections, these methods are briefly discussed from the perspective of system identification.

3.5.1 Minimum Manhattan Distance

If the identified Pareto front is bent, then the solutions in the ‘knee-region’ of the front represents the best overall compromise. Therefore, many existing a posteriori MCDM methods have been designed to select the knee-point solutions. Hence, in the absence of any preferences, the knee-point structures can be considered for further analysis. To this end, Minimum Manhattan Distance (MMD) approach has been selected, which is similar to Divide $\&$ Conquer approach to identify the knee-point solution, albeit MMD is computationally more efficient [67].

For the given set of non-dominated structures, $\Gamma$, MMD begins by identifying the hypothetical ideal point, $\vec{J}^{\star}$, in the corresponding objective space, as follows:

[TABLE]

Subsequently, for each non-dominated structure $\mathcal{X}_{i}\in\Gamma$, the Manhattan distance, $\mathcal{D}(\cdot)$, is evaluated with respect to $\vec{J}^{\star}$:

[TABLE]

Note that the Manhattan distance $\mathcal{D}(\cdot)$ is determined in the normalized objective space to avoid scaling issues. In MMD, as the name suggests, the structure with the minimum $\mathcal{D}(\cdot)$ is selected as the final choice. However, a slightly different approach is followed in this study, where each structure is ranked in the ascending order of $\mathcal{D}(\cdot)$, and the top five structures are selected for further analysis.

It is worth to note that MMD has explicitly been designed to discourage any preference specification from the DM. Further, the knee-point selection approach of MMD may not necessarily align with the preferences of the DM. Hence, MMD is not recommended when DM can specify preferences. For such scenario, Multi-criteria Tournament Decision (MTD) [66] may be more appropriate, which is discussed in the following.

3.5.2 Multi-criteria Tournament Decision

In practice, DM may have specific preferences which can somehow be embedded in the MCDM process. For instance, in system identification, the DM can assign higher priority to the number of terms, ‘$\xi$’, as a parsimonious description of the given system is usually preferred. Such preferences can easily be converted into a priority weight, ‘$w$’, which highlights the importance of a particular search objective as per DM’s preference. A detailed investigation for determining the priority weights from abstract DM preferences can be found in [65]. The priority weights for this study are determined following the guidelines in [65, 66], as discussed in Section 3.4.2.

The Multi-criteria Tournament Decision (MTD) approach essentially ranks the given non-dominated structures based on priority weights specified by the DM [66]. For this purpose, for each structure $\mathcal{X}_{i}\in\Gamma$, a pairwise comparison with the remaining structures is carried out, considering only the $p^{th}$ search objective, $p\in[1,m]$. This is accomplished by the following tournament function:

[TABLE]

Thus, for the $p^{th}$ objective, the tournament function $T_{p}(\mathcal{X}_{i},\Gamma)$ determines the total number of structures compared to which $\mathcal{X}_{i}$ yields better objective value. For the structure under consideration, the procedure is repeated to determine the tournament function across all the ‘$m$’ objectives. A similar procedure is repeated for each of the structure $\mathcal{X}_{i}\in\Gamma$. Subsequently, the following geometric aggregation function is used to derive a global rank which considers the tournament function for all the ‘$m$’ objectives as well as the priority weights which are specified by DM, $[w_{1},w_{2},\dots,w_{m}]$:

[TABLE]

It is easier to follow that a higher value of the global rank ‘$\mathcal{R}(\cdot)$’ is desirable. In this study, the structures are ranked in the descending order of $\mathcal{R}(\cdot)$. On the basis of this rankings, the top five structures are selected for further analysis.

4 Investigation Framework

In the following, the framework of this investigation is discussed. First, the search environment of this investigation is discussed in Section 4.1. The test non-linear systems being used to benchmark the search performance are discussed next in Section 4.2. Finally, the performance metrics being used to evaluate search from the quantitative and qualitative perspectives are discussed in Section 4.3 and 4.4.

4.1 Search Environment

The overall procedure followed for the multi-objective structure selection is outlined in Algorithm 3. The search process begins by formulating a priori preferences as goals. All the results in this study are obtained with the following goal specification: $\{\xi_{limit},\mathcal{E}_{limit}\}=\{20,30\}$. Subsequently, the non-dominated structures are gathered over multiple independent runs of the MOEA due to stochastic nature of the MOEAs. The dominance of the ensuing pool of the structures is again determined and only non-dominated structures are considered for further analysis. For the a posteriori decision making, the identified non-dominated structures are ranked using the Manhattan distance $D(\cdot)$ in MMD (see Section 3.5.1) and the rank $\mathcal{R}(\cdot)$ in MTD (see Section 3.5.2). The priority weights required for MTD approach are determined using preference ordering approach (see Section 3.4.2).

All of the compared MOEAs are implemented in MATLAB. For fair comparison, each MOEAs is terminated after 25,000 Function Evaluations (FEs). The population size, ‘$ps$’ is fixed to $50$ for all the MOEAs. The size of external archives in SPEA-II and MOEA/D is also set to $50$. In SPEA-II, the density estimate is determined with respect to $k^{th}$ nearest neighbor, where $k=10$. In MOEA/D, the sub-problems are created by uniformly spaced weight vectors in $[0,1]$. The neighborhood size, ‘$T$’, is fixed to $10\%$ of the population size. The new solution is allowed to replace a total of two solutions in the neighborhood, i.e., $nr=2$. To accommodate the stochastic nature of the MOEAs, the non-dominated structures over $40$ independent runs are selected.

Further, the uniform crossover is selected as the recombination mechanism in all the MOEAs on the basis of detailed comparative analysis (see Section 7.3). The following crossover and mutation probabilities for MOEAs are selected on the basis the performance sweet spot analysis discussed in Section 7.2: NSGA-II: $\{p_{c},\ p_{m}\}=\{0.9,0.006\}$; SPEA-II: $\{p_{c},\ p_{m}\}=\{0.7,0.008\}$; MOEA/D: $\{p_{c},\ p_{m}\}=\{0.8,0.008\}$.

4.2 Test Non-linear Systems

In this study, a total of eight known benchmark nonlinear systems have been included. These systems have been used by various researchers to study the efficacy of the structure selection approaches [13, 21, 68, 69, 36, 32]. The first six systems, shown in Table 1, are being used for comparative evaluation of MOEAs. The other known system has been included to highlight the exception to the specification of a priori goal point in Section 5.2.

For the comparative evaluation, a total of $1000$ data-points are collected for each system shown in Table 1. Each system is excited by a white noise input with either Gaussian (denoted by ‘WGN’) or uniform (denoted by ‘WUN’) distribution. The model set ($\mathcal{X}_{model}$) of $165$ NARX terms is obtained by the following specifications in NARX model (2): $[n_{u},n_{y},n_{l}]=[4,4,3]$. For structure selection purposes, $700$ data-points are used for the parameter estimation and the remaining data-points are used for the validation, i.e., $\mathcal{N}_{v}=300$.

In addition to these known systems, two case studies are also included: the identification of 1) Duffing’s oscillator and 2) nonlinear wave-forces acting on fixed structures. These case studies are discussed respectively in Section 5.3 and 5.4.

4.3 Quantitative Performance Metrics

The part of this study aims to compare the search performance of MOEAs on the identification problems, which requires a quantitative performance metric to infer the differences in the identified non-dominated structures. Further, a quantitative metric is essential to investigate the effects of search parameters and the choice of the crossover operator which will be discussed in Section 7.1. Usually, MOEAs are evaluated on the basis of the following two criteria: 1) proximity of the identified solutions to the optimal Pareto front and 2) diversity of the identified solutions. Most of the existing performance metrics have been designed to integrate the evaluation of these criteria to yield a singular value [54, 70]. In this study, two of such performance metrics are being used to evaluate the performance of MOEAs: Set Coverage Metric and Hyper-volume Indicator [71], which are discussed briefly in the following:

Consider a structure selection problem with ‘$m$’ number of objectives being solved by two MOEAs, ‘$A$’ and ‘$B$’. Let the set of non-dominated structures and the corresponding set of objective function vectors identified by these algorithms be denoted by $\Gamma_{A}$, $\Gamma_{B}$, $\Lambda_{A}$ and $\Lambda_{B}$, i.e.,

[TABLE]

The set coverage metric, $\mathcal{C}(\cdotp)$, is designed to compare two sets of non-dominated structures. For instance, for $\Gamma_{A}$ and $\Gamma_{B}$, it essentially determines the number of structures in $\Gamma_{B}$ which are dominated by the structures in $\Gamma_{A}$, as follows [71]:

[TABLE]

It is easier to follow that this metric is bounded between $[0,1]$. It is worth emphasizing that the $\mathcal{C}$-metric may not be symmetric [71]. Hence, the comparison of $\Gamma_{A}$ and $\Gamma_{B}$ requires evaluation of both $\mathcal{C}(A,B)$ and $\mathcal{C}(B,A)$. It is obvious that $\mathcal{C}(A,B)>\mathcal{C}(B,A)$ implies the search performance of algorithm ‘$A$’ is better than that of Algorithm ‘$B$’, i.e., $\Gamma_{A}\preceq\Gamma_{B}$.

In contrast to $\mathcal{C}$-metric, the hyper-volume (HV) indicator is designed to independently evaluate the given set of non-dominated structures. In essence, the HV-indicator measures the region in the objective space which is dominated by the given set of non-dominated structures, say $\Gamma_{A}$, and bounded by a reference point, $\vec{r}\in\mathbb{R}^{m}$, which is given by,

[TABLE]

A higher value of this metric denotes a better performance. In this study, HV-indicator is evaluated following the dimension sweep algorithm [72, 73]. Further, to avoid scaling issues, it is determined using the normalized objective values. Note that HV-indicator may be sensitive to the specification of the reference point, $\vec{r}$. To avoid such issues, we follow the guidelines in [74] to determine $\vec{r}$.

4.4 Qualitative Evaluation of Search Outcomes

It is worth emphasizing that a discrete NARX model of a continuous-time system is not unique, i.e., there may exist several valid NARX structures for the system under consideration [4]. This issue is further explored later through the case studies in Section 5.3 and 5.4. However, this is not the case for the known discrete non-linear system shown in Table 1, which are being used for the comparative evaluation. Given that the identification data are directly generated from these systems, there exists a unique and known structure for each system in Table 1. Hence, for these systems, the non-dominated structures identified by MOEAs can easily be evaluated. To this end, each non-dominated structure is ‘qualitatively’ evaluated following the search outcome definitions developed in our previous investigations [25, 53].

Consider the following subsets of NARX terms for this purpose,

• •

$\varnothing$ : the null set

• •

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

• •

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

• •

$\mathcal{X}_{k}$ : $k^{th}$ non-dominated structure identified by an MOEA, $\mathcal{X}_{k}\in\Gamma^{\ast}$

• •

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

Following these definitions, each non-dominated structure can be categorized into one of the following search outcome [25]:

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., $\mathcal{X}_{k}=\mathcal{X}^{\star}$ and $\mathcal{X}_{k}^{spur}=\varnothing$ 2. 2.

Over Fitting :

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

Under Fitting-1 :

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

Under Fitting-2 :

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

5 Results: Identification Perspective

This study aims to establish the multi-objective evolutionary algorithms (MOEAs) as viable alternatives for structure selection. To this end, several issues encountered in multi-objective structure selection are investigated. First, the identification of a set of benchmark nonlinear systems is discussed in Section 5.1. Next, the issues related to the specification of goal point are discussed through an illustrative example in Section 5.2. Finally, two case studies are included in Section 5.3 and 5.4. The identification of discrete-time model of the continuous-time Duffing’s oscillator is taken as the first case study. Next, the practical problem of non-linear wave-force identification is considered. The objective of these case studies is to highlight key distinctive ability of MOEAs, i.e., identification of multiple valid system structures.

5.1 Identification of Discrete-time Systems

The first phase of the investigation begins by considering several discrete-time nonlinear systems shown in Table 1. The non-dominated structures for these systems are identified using all the MOEAs and as per the search environment described in Section 4.1 and Algorithm 3. It is worth emphasizing that, since the structure of these systems are known, the correctness of the proposed approach can be easily validated. As discussed and defined in Section 4.4, the non-dominated structures identified by MOEAs are categorized using qualitative search outcomes. Since the MOEAs are designed to obtain diverse solutions across the entire Pareto front, the approximate Pareto set may contain several under-fitted and/or over-fitted structures.

From the system identification perspective, the ideal search outcome is ‘exact-fitting’. However, the over-fitted structures can also be tolerated provided the spurious terms can be removed through some refinement procedure. To investigate this further, the identified non-dominated structures are refined through a simple null-hypothesis test [18] on the coefficients. Subsequently, the refined structures are categorized into qualitative search outcomes as per the definitions in Section 4.4. Following these steps, the frequency of each search outcome is determined for the approximated Pareto set, as shown in Table 2. For instance, the approximate Pareto set identified by NSGA-II for the system $\mathcal{S}_{2}$ contains a total of $17$ non-dominated structures. Each of these structures is refined using null-hypothesis test and categorized into one of the four possible search outcomes. Following these steps, $16$ non-dominated structures are categorized as ‘exact-fitting’ whereas only one non-dominated structure is classified as ‘under-fitting’ (as seen in Table 2).

It can be argued that ‘success’ of the given MOEA can be determined on the basis of total number of ‘exact-fitting’ outcomes. Pursuing this argument, it is interesting to see that most of the identified non-dominated structures are categorized as exact-fitting (see Table 2), which implies that these structures could include all the system terms. It is, therefore, safe to conclude that all the MOEAs have been successful in identifying the correct structure for all the benchmark discrete-time systems in Table 1.

5.2 Comments on the a priori Goal Point

This study takes a pragmatical approach to set the a priori goal point, i.e., $\{\xi_{limit},\mathcal{E}_{limit}\}=\{20,30\}$, as discussed in Section 3.4.1. In particular, these specifications have been selected to guide the search towards structures with fewer terms. However, it is important to note that these specifications are ultimately system-dependent. In the following, a numerical example is considered to highlight this issue along with a simple approach to handle such exceptions.

Consider the following non-linear system with $23$ number of terms:

[TABLE]

A total of $1000$ input-output data points are collected by exciting the system with $u\sim WUN(0,1)$ and $e\sim WGN(0,0.01^{2})$. The model set ($\mathcal{X}_{model}$) of $286$ NARX terms is obtained with the following specification in (2): $[n_{u},n_{y},n_{l}]=[5,5,3]$.

NSGA-II is applied to identify this system with $\{\xi_{limit},\mathcal{E}_{limit}\}=\{20,30\}$, as per the procedure outlined in Section 4.1 and Algorithm 3. The approximate Pareto front (APF) for this scenario is shown in Fig. 3(a). Note that the specified cardinality limit in this case is less than the actual system cardinality, i.e., $\xi_{limit}<\xi^{\ast}$. Hence, the search is directed towards only the under-fitted structures (see Section 4.4). This is not apparent in APF shown in Fig. 3(a). However, this can easily be detected by observing the behavior of classical Information Criteria (IC) (e.g., BIC, LILC) for the non-dominated structures. For instance, Fig. 3(b) shows the IC obtained with the non-dominated structures. It is clear that the minimum value for both BIC and LILC is obtained at $\xi=20$. However, it is not clear whether this is the actual knee-point/plateau of these criteria. Hence, the cardinality preference is updated to $\xi_{limit}=50$ and the identification procedure is repeated. The consequent Pareto front and IC are shown respectively in Fig. 3(c) and 3(d). With the updated preference, a continual increase in information criteria is observed for $\xi\geq 24$. This, therefore, confirms that the plateau observed in the information criteria for $\xi\in[23,27]$ indeed points towards the true system cardinality. Hence, no further increase in $\xi_{limit}$ is required.

As illustrated by this example, evaluation of the specified a priori preferences and any required updates can easily be accomplished by observing the behavior of the information criteria of the non-dominated structures.

5.3 Discrete Models for Nonlinear Continuous-Time Systems

In this part of the investigation, identification of a discrete-time model for a continuous-time system is considered as a case study to evaluate the multi-objective framework further. In particular, the NARX model is fitted to identify the dynamics of Duffing’s oscillator, which is given by,

[TABLE]

For this purpose, a total of $1000$ data points are generated by exciting the system with a uniform white noise sequence, i.e., $u(\cdotp)\sim WUN(0,1)$, and sampling the output at $500\ Hz$. The system parameters are set as follows: $\omega_{n}=45\pi$, $\zeta=0.01$ and $\varepsilon=3$. The model set ($\mathcal{X}_{model}$) with a total of $286$ terms is obtained by specifying $[n_{u},\ n_{y},\ n_{l}]=[5,\ 5,\ 3]$ in (2). A total of $700$ data points is used for estimation purposes and the remainder is used for validation.

The NARX models for the Duffing’s oscillator are identified following the procedure outlined in Section 4.1 and Algorithm 3. Note that, in the following, only the results obtained by SPEA-II are discussed, as it yields comparatively better approximation of the Pareto set (to be discussed in Section 6). A total of $17$ non-dominated structures are identified by SPEA-II. Subsequently, the identified non-dominated structures are ranked based on Manhattan distance ‘$\mathcal{D}(\cdot)$’ (see Section 3.5.1) and the global rank ‘$\mathcal{R}(\cdot)$’ (see Section 3.5.2). On the basis of these rankings, three non-dominated structures are selected for further analysis. Table 3 shows the objective functions of the selected non-dominated structures along with the corresponding Manhattan distance and the global rank. Note that predictive performance, shown in Table 3, is determined over the model-predictive output. The selected structures and the corresponding coefficients are given by the following models:

[TABLE]

It is clear that all the identified models could yield practically acceptable predictive performance and could be used to describe the dynamics of the Duffing’s oscillator. This is expected as the discrete model for the continuous-time system is not unique [4].

To investigate whether the identified models indeed describe the same dynamics, the Generalized Frequency Response Functions (GFRFs) are computed from these models by the harmonic probing procedure [12] and are compared with that of the continuous-time model. Note that the GFRFs are invariant descriptors of the system dynamics. This implies that if the identified models are correct the corresponding frequency response function should match with that of the continuous-time model.The GFRFs can easily be determined for differential equations models, the polynomial NARX models and the other system representations by mapping them into the frequency domain [4, 12, 75].

The linear and the third order GFRFs computed from the continuous-time model and the identified discrete models are shown in Fig. 4-6. We begin by identifying the unique frequency domain ‘traits’ of the Duffing’s oscillator from the GFRFs shown in Fig. 4(a), 5(a) and 6(a). As seen Fig. 4(a), the linear resonant frequency of Duffing’s oscillator is at $22.5\ Hz$. Further, the third order GFRF, in Fig. 5(a) and 6(a), is characterized by the presence of several peaks and ridges. The peaks occur whenever the system is excited by a three tone signal with $f_{1}=f_{2}=f_{3}=\pm 22.5\ Hz$. Further, the ridges occur in either of the following conditions: $f_{1}$, $f_{2}$ or $f_{3}=\pm 22.5\ Hz$ or $f_{1}+f_{2}+f_{3}=\pm 22.5$.

It is interesting to see that all the identified models ($\mathcal{M}_{D_{1}}$, $\mathcal{M}_{D_{2}}$, $\mathcal{M}_{D_{3}}$) could replicate the aforementioned frequency domain behavior of the Duffing’s oscillator (as seen in Fig. 4-6). Hence, this gives conclusive evidence that these apparently different models indeed describe the dynamics of the Duffing’s oscillator. However, since the model $\mathcal{M}_{D_{1}}$ enables relatively compact system description, it is recommended following the principle of parsimony.

5.4 Identification of nonlinear wave forces

Next, we focus on a practical problem from the field of computational hydrodynamics which involves modelling of wave forces acting on fixed bodies in a viscous fluid. In particular, the objective is to capture the wave force dynamics from the flow velocity (input, $u$) and the wave force (output, $y$) histories [76]. A total of $1400$ data points are used for the identification purposes, which are shown in Fig. 7. Further details about this dataset can be found in [76].

The models for wave-force dynamics are identified following the procedure outlined in Section 4.1 and Algorithm 3. To this end, the model set ($\mathcal{X}_{model}$) of $286$ NARX terms is generated by the following specification in (2): $[n_{u},n_{y},n_{l}]=[5,5,3]$. A total of $1000$ data points are being used for estimation purposes and the rest are used for validation. With these specifications, the approximate Pareto fronts (APF) are shown in Fig. 7(c). It is clear that NSGA-II could yield better structures for the wave-force data (to be discussed in Section 6). Hence, only the results obtained by NSGA-II are discussed in the following.

A total of $10$ non-dominated structures are identified by NSGA-II. These structures are ranked based on Manhattan distance ‘$\mathcal{D}(\cdot)$’ (see Section 3.5.1) and the global rank ‘$\mathcal{R}(\cdot)$’ (see Section 3.5.2). Based on these rankings, following three models are selected for further analysis:

[TABLE]

Table 4 gives the objective functions of the selected non-dominated structures along with the corresponding Manhattan distance and the global rank. The prediction errors (NMSE and MSE) in Table 4 are determined over the model predicted output with the validation data-points, which are shown in Fig. 8. It is clear that there is no significant difference in the predictive performance of the identified models.

For further validation, the frequency domain behavior of the identified models is compared using GFRFs. Note that even though the system dynamics may be explained by distinct discrete-time models, they are uniquely defined in the frequency domain [4]. Hence, if the apparently different models have captured the wave force dynamics then the corresponding GFRF should match. To investigate this, GFRFs of the identified models are determined by mapping them to frequency domain, and are shown in Fig. 9. These results do not show any significant difference in the linear resonant frequency or the peak magnitude.

Further, it is interesting to note that the similar frequency response behavior was observed in the distinct NARX model identified for the same wave force data in our previous investigation [25].

6 Results: Search Perspective

The objective of this part of the study is to determine whether there exists any significant difference in the search performance of the MOEAs. To this end, in the following, MOEAs are first compared using the quantitative performance metrics, drawing on the recommendations in [70]. Subsequently, a rigorous statistical analysis is carried out to determine the significance of the results.

For this test, in addition to the seven discrete-time systems ($\mathcal{S}_{1}$-$\mathcal{S}_{7}$), a nonlinear continuous-time system (Duffing’s oscillator, Section 5.3) and a practical case study (wave-force identification, Section 5.4) are included. For each test system, the non-dominated structures are identified over $40$ independent runs of MOEAs as per the search environment described in Section 4.1 and Algorithm 3. The goal point for the system $\mathcal{S}_{7}$ is updated to $\{\xi_{lim},\ \mathcal{E}_{lim}\}=\{50,\ 30\}$, as discussed in Section 5.2. The set of non-dominated structures which are obtained over $40$ runs is treated as the approximate Pareto set (APS). Further, for a particular system, the APS identified by NSGA-II, SPEA-II and MOEA/D are respectively denoted by $\Gamma_{A}$, $\Gamma_{B}$ and $\Gamma_{C}$.

6.1 Dominance based comparison of MOEAs

First, the set coverage metric $\mathcal{C}(\cdot)$ is considered which allows the pairwise comparison of the given APS. Note that, since $\mathcal{C}(\cdot)$ is not symmetric, the comparison of given APS (say, $\Gamma_{A}$ and $\Gamma_{B}$) requires evaluation of both $\mathcal{C}(\Gamma_{A},\Gamma_{B})$ and $\mathcal{C}(\Gamma_{B},\Gamma_{A})$. Hence, for sake of clarity, the difference in these two metrics $\Delta\mathcal{C}(\cdot)$ is reported in Table 5. This is defined as:

[TABLE]

A positive value of $\Delta\mathcal{C}(\Gamma_{A},\Gamma_{B})$ implies that $\Gamma_{A}$ weakly dominates $\Gamma_{B}$ whereas a negative value implies otherwise. This metric is determined for all the test systems and it is shown in Table 5. The results are shown for three pair-wise comparisons among NSGA-II, SPEA-II and MOEA/D. Based on these pair-wise comparisons, for each system, the dominance ordering for MOEAs is determined and is shown in the last column of Table 5. For instance, for the system $\mathcal{S}_{1}$, when NSGA-II is compared to SPEA-II and MOEA/D, a positive value of $\Delta\mathcal{C}(\cdot)$ is obtained, i.e., $\Gamma_{A}\preceq\Gamma_{B}$ and $\Gamma_{A}\preceq\Gamma_{C}$. Further, for the same system, when SPEA-II is compared to MOEA/D, a negative value of $\Delta\mathcal{C}(\cdot)$ is obtained, i.e., $\Gamma_{C}\preceq\Gamma_{B}$. Hence, based on these pairwise comparisons, the following dominance ordering of MOEAs can be inferred for $\mathcal{S}_{1}$ : $\Gamma_{A}\preceq\Gamma_{C}\preceq\Gamma_{B}$.

By inspecting the dominance relations shown in Table 5, a clear pattern emerges. For instance, NSGA-II outperforms both SPEA-II and MOEA/D on most of the test systems (six out of nine). Further, SPEA-II exhibits better performance than MOEA/D on six test systems. Hence, considering all the test systems, the performance of MOEAs can be ranked as: NSGA-II$>$SPEA-II$>$MOEA/D.

6.2 Non-parametric statistical comparison of MOEAs

Owing to the stochastic nature of MOEAs, further statistical analysis is carried out to determine the significance of the comparative analysis. In particular, the objective is to test the hypothesis that NSGA-II is significantly better than the other two MOEAs, SPEA-II and MOEA/D. For this purpose, non-parametric statistical tests are considered; as in practice the results obtained by MOEAs usually do not satisfy the necessary conditions of parametric tests [77]. Therefore, in the following, multiple non-parametric statistical comparisons are carried out based on the recommendations in [78, 77].

Note that while the coverage metric $\mathcal{C}(\cdot)$ determines the dominance ordering of the given approximate Pareto set (APS), it does not provide any information on the degree of improvement. Hence, for the statistical analysis, the hyper-volume (HV) metric is considered (see Section 4.3). The reference point necessary for HV calculation is determined following the guidelines in [74]. Prior to the determination of HV, the objective values are normalized and bounded in $[0,1]$. Following this procedure, HV is determined for all the test systems and are shown in Table 6. The significance of these results is determined in two steps.

First, the Friedman Two-way Analysis of Variance by Ranks [78, 77] is applied to test the null hypothesis that median value of HV obtained with MOEAs are equal. For this purpose, the hyper-volume values are ranked from ‘$1$’ (highest) to ‘$3$’ (lowest). The procedure is repeated for all the test systems. Subsequently, the average value of ranks, the test statistic and the corresponding p-value are shown in Table 6. These results indicate that the null-hypothesis can easily be rejected for $95\%$ confidence interval. This implies that there exists a significant difference in the performance of MOEAs.

Next, multiple comparisons of NSGA-II with the other MOEAs is carried out through a set of interconnected hypotheses. In particular, the determination of two interconnected null hypotheses is required to compare NSGA-II with SPEA-II and MOEA/D. Each null hypothesis ($H_{0}$) implies that the compared MOEA is significantly better than NSGA-II. The test statistic and $p-$value needed to evaluate these hypotheses are determined following the guidelines in [77, 78]. The $p-$value are adjusted for multiple comparisons following the Hommel’s post-hoc procedure [77], and denoted as Adjusted p-values (APV). The hypotheses are evaluated at $95\%$ confidence interval. The outcomes of this test are shown in Table 7. This test disproves both the null-hypotheses and thereby confirms that NSGA-II is significantly better than SPEA-II and MOEA/D.

7 Robustness to Control Parameters

A judicious tuning of search parameters is crucial to the performance of any evolutionary algorithm (EA) [79]. For instance, the cross-over ($p_{c}$) and mutation probabilities ($p_{m}$), the choice of recombination mechanism such as single-point or uniform crossover can have significant effect on the search performance. These parameters can be categorized into quantitative or qualitative parameters.

Although there exist several parameter tuning approaches, e.g. see [79], in most of the existing investigations, the parameters are selected either based on the prevalent practice or based on limited trial-and-error approach [41, 42, 43, 46, 45, 49, 47, 48]. Even though MOEAs considered in this study have been developed over past two decades, their search performance has received relatively less attention from the identification perspective. It is therefore crucial to evaluate the robustness of these algorithms with respect to quantitative and qualitative search parameters. This has been the motivation behind the empirical investigation which is discussed in the following.

The quantitative control parameters include but not limited to the size of the population and external archive (if any), the cross-over ($p_{c}$) and mutation probabilities ($p_{m}$), and other algorithm specific parameters, e.g., neighborhood size ‘$T$’ in MOEA/D. In this study, the effects of the cross-over ($p_{c}$) and mutation probabilities ($p_{m}$) on the search performance are of particular interest. For this purpose, the search performance of MOEAs is determined over a grid in the parameter space. Each point on this grid denotes a particular combination of parameters $\{p_{c},\ p_{m}\}$. For each parameter combination, the search performance is quantified using the hyper-volume metric. This will determine the region of better performance over the grid-map. Such promising regions in the parameter space are often referred to as ‘sweet-spot’ [50, 51]. It is easier to follow that the robustness of the given MOEA, with respect to the search parameter, can indirectly be estimated from the sweet-spot in the parameter space. More details and the results of this experiment are respectively discussed in Section 7.1 and 7.2.

The qualitative parameters mainly reflect the choice of algorithm design, e.g., choice of selection, recombination and mutation mechanisms. In MOEAs, the proximity and diversity preserving mechanisms are closely tied to the selection process. The choice of selection mechanism is thus dependent on the framework of MOEA, e.g., crowded tournament selection in NSGA-II. Hence, for each MOEA, the selection mechanism is chosen in the original research is considered. In contrast, the choice of the recombination and the mutation mechanisms is not dependent on the framework of MOEA. Hence, it is possible to incorporate suitable mechanism. In particular, we are interested in evaluating the effects of recombination mechanism, which will be discussed in Section 7.3.

7.1 Experimental Setup

The search setup for the robustness analysis is outlined in Algorithm 4. For this test, the grid is designed in the parameter space using $10$ distinct values of $p_{c}$ and $p_{m}$, where $p_{c}\in[0.1,1]$ and $p_{m}\in[0.001,0.02)$. This gives a total of $100$ distinct combinations of $\{p_{c},p_{m}\}$. The search performance is evaluated on all the discrete time systems ($\mathcal{S}_{1}-\mathcal{S}_{7}$). In Algorithm 4, the total number of parameter combinations and systems are respectively denoted by ‘$n_{comb}$’ and ‘$n_{sys}$’, i.e., $n_{comb}=100$ and $n_{sys}=7$.

The goal of this test is to determine the unary performance metric corresponding to each parameter combination on the test grid. To this end, the test begins by selecting the MOEA and the crossover mechanism, and it is further carried out in the following two stages: 1) Search for Approximate Pareto Set (APS) and 2) Evaluation of Performance Metric.

In the first stage, for each test system, APS is identified corresponding to each parameter combination. The steps involved in this stage are outlined in Line 4-4, Algorithm 4. Due to stochastic nature of MOEAs, the non-dominated structures identified over ‘$R$’ independent runs ($R=20$) are determined to be APS. For instance, in Algorithm 4, ‘$\Gamma_{i,j}$’ and ‘$\Lambda_{i,j}$’ respectively denote APS and the corresponding front identified for the $i^{th}$ system with the $j^{th}$ parameter combination in MOEA. This procedure is repeated for each system across all the combination of parameters.

In the second stage, the identified APS are quantified using unary performance metric. The overall procedure involved in this stage is outlined in Line 4-4, Algorithm 4. In particular, for the given combination of parameters, the goal is to determine an aggregated value of performance metric across all the test systems. For this purpose, the hyper-volume ratio (HVR) is appropriate as it yields normalized values which are bounded in $[0,1]$. As the name suggests, HVR is in essence ratio of hyper-volume values of the identified APF and the ‘ideal’ Pareto front. Note that since the ‘ideal’ Pareto front is not known, it is approximated from the union set of APF which are identified with different parameter combinations, as outlined in Line 4-4, Algorithm 4. Subsequently, HVR is determined for each system across all the parameter combination from the corresponding APS. For example, in Algorithm 4, ‘$PM_{i,j}$’ denotes the HVR for the $i^{th}$ system with the $j^{th}$ parameter combination. Finally, for each parameter combination, the aggregated performance metric ($\overline{PM}$) is determined as shown in Line 4, Algorithm 4.

The overall procedure outlined in Algorithm 4 is repeated for all combinations of MOEAs (NSGA-II, SPEA-II and MOEA/D) and recombination mechanisms (Single Point and Uniform).

7.2 Sweet-spots in the Parameter Space

The objective of this part of the study is to determine the robustness of MOEAs to parameter specification. As discussed earlier, this can indirectly be determined by observing the region of ‘good’ performance in the search parameter space, i.e., the bigger the region of high performance, the higher is the robustness. To this end, an aggregated performance metric ($\overline{PM}$) across all the test systems is determined for each parameter combination, as per the procedure outlined in Algorithm 4. Following this procedure a contour map of performance($\overline{PM}$) is determined in the parameter space; as each parameter combination is essentially a point on the parameter grid (as discussed in the previous section). Consequently, the contours are used to highlight the region of high performance, which are often referred to as ‘sweet-spots’. Note that this approach is somewhat similar to earlier investigation of Goldberg [50] and Purshouse and Fleming [51], where the objective was to determine the performance sweet-spots over the grid of some quantitative or qualitative search parameters.

The results of this investigation are shown in Fig. 10. For each contour-map the corresponding color-bar indicates the scale of $\overline{PM}$. A higher value of $\overline{PM}$ is desirable which is represented by a lighter shade of yellow. Hence, the sweet-spot is highlighted by a yellow bounded area in each contour-map. It is interesting to see that the sweet-spots of NSGA-II and MOEA/D inflate/grow when uniform crossover is used (see Fig 10(a)-10(b), NSGA-II and Fig. 10(e)-10(f), MOEA/D), which implies that the robustness of NSGA-II and MOEA/D is better with uniform crossover (to be discussed in the next section). Hence, we focus only on the results obtained with the uniform crossover, i.e., Fig. 10(b) (NSGA-II), Fig.10(d) (SPEA-II) and Fig.10(f) (MOEA/D).

It is clear that the variance in $\overline{PM}$ with NSGA-II and SPEA-II is relatively low. Hence, the search performance will not significantly deteriorate even with improper choice of parameters. In contrast, the parameter selection is a critical issue in MOEA/D which is indicated by a high variance in $\overline{PM}$. From the contours, the sweet-spots of the MOEAs are determined as follows: $p_{c}\in(0.2,1]$ and $p_{m}\in[0.002,0.009]$ (NSGA-II, Fig. 10(b)); $p_{c}\in[0.3,1]$ and $p_{m}\in[0.003,0.01]$ (SPEA-II, Fig.10(d)) and $p_{c}\in[0.5,1]$ and $p_{m}\in[0.003,0.013]$ (MOEA/D, Fig.10(f)). Note that these performance sweet-spots can also be used as rule of thumb or a starting point to select the parameters for the other identification problems.

The comparative analysis suggests that the search performance of MOEA/D is poorer compared to the other dominance based MOEAs such as NSGA-II and SPEA-II. This could be ascribed to the choice of the control parameters. In this study, a detailed investigation is carried out to select the crossover and mutation probabilities as well as the recombination mechanism. However, further investigation is required to determine the sensitivity of MOEA/D to the other control parameters such as the choice of aggregation function, neighborhood size, ‘$T$’, and replacement neighbor size,‘$nr$’.

7.3 Recombination: Single Point or Uniform?

The objective of this part of the study is to determine if there exist any significant effects of the recombination mechanism on the search performance of MOEAs.

The motivation for this test is in part based on the previous investigation by the authors in which Simple Genetic Algorithm (SGA) was found to be more susceptible to the selection of spurious terms [52]. While the cause of this poor performance is still to be thoroughly investigated, in theory, it can be ascribed to poor selection and/or recombination mechanism. Especially, the single point crossover in SGA is designed to preserve low-order schemas which are considered to be the building blocks of the string [56]. This assumption, however, may not be valid in the structure selection, since the significant/system terms may not be positioned close to each other in the binary string representation, i.e., the desired schemas may not necessarily be of low-order. Further, it has been argued that the pursuit to preserve low-order schemas may lead to selection of bits that are not part of the desired schema [80]. This could explain the selection of spurious terms (bits) in our previous investigations. It has been argued that the parameterized uniform crossover can overcome these issues [80, 81]. Since MOEAs being used are in essence based on genetic algorithms, it is pertinent to evaluate the effects of recombination mechanisms on the search performance of MOEAs.

In particular, the goal here is to test the hypothesis that the uniform crossover can significantly improve the search performance of MOEAs in comparison to the single point crossover. To this end, the test outlined in Algorithm 4 is carried out with all the MOEAs and in conjunction with both Single Point and Uniform crossovers. Following these steps, the performance metric ($PM$) is obtained for all possible instances, where each instance denotes a unique grouping of MOEA, recombination mechanism and parameter combination. For sake of clarity, the results with three parameter combinations (out of 100) for $\mathcal{S}_{1}$ are shown in Table 8. These results can be interpreted as follows: Consider the first row in Table 8. These results indicate that when the crossover and mutation probabilities are set to $\{0.1,0.001\}$, NSGA-II can yield better results with uniform crossover; since the corresponding $PM$ value ($0.9764$, Table 8) is higher than that of the single point crossover (i.e., $0.9696$, Table 8). Following the same arguments, local conclusions can be inferred for other parameter combinations. It can be argued that the crossover mechanism, which yields better performance across higher number of parameter combination, is desirable. Hence, for each MOEA, the following hypotheses are evaluated using the Wilcoxon matched-pairs signed-ranks test [78, 77]:

• •

Null Hypothesis ($H_{0}$): The median value of performance metrics from single point and uniform crossover are equal, i.e., $H_{0}\mathrel{\mathop{\mathchar 58\relax}}\mu_{D}=0$

• •

Alternative Hypothesis ($H_{1}$): The median value of performance metric from uniform crossover is higher than that of the single point crossover, i.e., $H_{1}\mathrel{\mathop{\mathchar 58\relax}}\mu_{D}>0$

These hypotheses are determined for each test system using the performance metric ($PM$) values across all the parameter combinations. The outcome of the hypotheses evaluation at $95\%$ confidence interval are shown in Table 9. It is interesting to see that, with MOEA/D, the null or equivalence hypothesis is rejected for all the test systems, i.e., there is enough evidence to support the argument that MOEA/D would yield better search performance with the uniform crossover. In contrast, the null hypothesis is rejected for few systems with the other MOEAs (4 with NSGA-II and 3 with SPEA-II). This indicates that with NSGA-II and SPEA-II, uniform crossover will provide better or equivalent performance to single point crossover. Hence, on the basis of these observations, the uniform crossover is selected in this study.

8 Conclusions

A Multi-Objective Structure Selection (MOSS) framework for the identification of nonlinear system is proposed. In particular, the structure selection is formulated as a Multi-Criteria Decision Making (MCDM) problem to balance the bias-variance dilemma. The proposed framework is critically analyzed both from the perspective of system identification and search algorithm through comparative analysis on several benchmark nonlinear systems and a practical identification problem. A thorough parameter sensitivity analysis has also been carried out to determine the robustness of MOEAs. These results show that the region of high performance or sweet spots in the parameter space is very large, especially for the dominance based MOEAs. This implies that the choice of the control parameters are not critical to the search performance of MOEAs, provided these lie in the identified sweet spots.

The key benefit of MOSS is the ability to identify multiple valid discrete-time models for continuous-time systems. This has been demonstrated by validating the identified models in frequency domain using generalized frequency response functions.

The proposed MOSS framework can easily accommodate distinct MOEA architectures. In this investigation, this is demonstrated by dominance based and decomposition based architectures. In addition, it is possible to incorporate distinct preferences to focus the search in the region of interest to the decision maker. The proposed approach therefore provides a flexible generic identification framework.

Acknowledgement

Faizal Hafiz is grateful to Education New Zealand for supporting this research through New Zealand International Doctoral Research Scholarship (NZIDRS).