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.


Making Agents' Abilities Explicit

Alternating-time temporal logics (ATL/ATL*) represent a family of modal logics for reasoning about agents' strategic abilities in multiagent systems (MAS). The interpretations of ATL/ATL* over the semantic model Concurrent Game Structures (CGS) usually vary depending on the agents' abilities, for instance, perfect vs. imperfect information, perfect vs. imperfect recall, resulting in a variety of variants which have been studied extensively in literature. However, they are defined at the semantic level, which may limit modeling flexibilities and may give counter-intuitive interpretations. To mitigate these issues, in this work, we propose to extend CGS with agents' abilities and study the new semantics of ATL/ATL* under this model. We give PSACE/2EXPTIME model-checking algorithms for ATL/ATL* and implement them as a prototype tool. Experiment results show the practical feasibility of the approach.

Introduction
Multiagent systems (MAS) consisting of multiple autonomous agents are a wide adopted paradigm of intelligent systems. Game-based models and associated logics, as the foundation of MAS, have received tremendous attention in recent years. The seminar work [Alur et al., 2002] proposed concurrent game structures (CGS) as the model of MAS and alternating-time temporal logics (typically ATL and ATL * ) as specification languages for expressing temporal goals. In a nutshell, a CGS consists of multiple players which are used to represent autonomous agents, components and the environment. The model describes how the MAS evolves according to the collective behavior of agents. ATL/ATL * , an extension of the Computational Tree Logics (CTL/CTL * ), features coalition modality A . The formula A ϕ expresses the property that the coalition (i.e. the agent group) A has a collective strategy to achieve a certain goal specified by ϕ.
A series of extensions of ATL-like logics have been studied which take different agents' abilities into account. These abilities typically include whether the agents can identify the current state of the system completely or only partially (perfect vs. imperfect information), and whether the agents can memorize the whole history of observations or simply part of them (perfect vs. imperfect recall). Different abilities usually induce distinct semantics of logics, which are indeed necessary because of the versatility of problem domains. These semantic variants and their modelchecking problems comprise subjects of active research for almost two decades, to cite a few [Schobbens, 2004;Jamroga and van der Hoek, 2004;Ågotnes et al., 2007;Dima and Tiplea, 2011;Bulling and Jamroga, 2011;Laroussinie and Markey, 2015].
While the agents' abilities play a prominent role [Bulling and Jamroga, 2014], the semantics of ATL-like logics only refers to them implicitly. In other words, the logic per se does not specify what ability an agent has; instead one can infer the ability an agent requires by examining the goal specified in the logic. This approach, being elegant and valuable to understand the relationship between different abilities, suffers from a few shortcomings: (1) From the modelling perspective, it is common in practice that agents in an MAS vary in their abilities (for instance, agents modeling sensors may not identify the complete state of system so can only use strategies with imperfect information). When building up a model, these abilities ought to be encoded explicitly. Such modeling flexibility is not supported by existing formalisms. (2) From the semantic perspective, ATL-like logics may exhibit some counter-intuitive semantics. A ϕ is interpreted as the coalition A has a collective strategy to achieve the goal ϕ "no matter what the other agents do" rather than "no matter which strategies the other agents choose", hence, neglects the (multi-player) game nature in the evolution of MAS. For instance, when it comes to the imperfect information/recall setting, only agents in the formula A ϕ are assumed to use imperfect information/recall strategies, while the agents not in A may still use perfect information and perfect recall strategies. Even worse, if the coalition modalities are nested, the same agent may have different abilities to fulfill the objectives specified in different sub-formulae, which results in inconsistency. This phenomenon has also been mentioned in [Mogavero et al., 2014;Cermák et al., 2014] which proposed a strategic logic does make explicit references to strategies of all agents including those not in A. However, all agents in strategic logic should have same abilities.
To summarize, it occurs to us that the current approach in which the temporal formulae are with implicit agents' abilities at the semantics level impedes necessary modeling flexibility and often yields unpleasant (even weird) semantics. Instead, we argue that coupling agents' abilities at the syntactic level of system models would deliver a potentially better approach to overcome the aforementioned limitations. Bearing the rationale in mind, we propose a new MAS model, named Agents' Abilities Augmented Concurrent Game Structures (ACGS), which encompasses agents' abilities explicitly. We investigate ATL/ATL * over ACGS. We show that in general the new semantics of ATL/ATL * over ACGS is incomparable with others even if the underlying CGS models are the same. We also study the model-checking problem of ATL/ATL * over ACGS. We show that this problem is generally undecidable, as the problem of ATL over CGS under the imperfect information and perfect recall setting is already undecidable [Dima and Tiplea, 2011]. However, we manage to show that the model-checking problem for ATL * (resp. ATL) on ACGS is 2EXPTIME-complete (resp. in PSPACE) when the imperfect information and perfect recall strategies are disallowed. We implement our algorithms in a prototype tool and conduct experiments on some standard applications from the literature. The results confirm the feasibility of our approach.
The source code of our tool is available at [MCMAS-ACGS, 2018] which also includes some further experiments and comparison of ATL/ATL * semantics between CGS and ACGS.

Concurrent Game Structures
Given an infinite word ρ = s 0 s 1 . . ., we denote by ρ j the symbol s j , ρ [0.. j] the prefix s 0 . . . s j , and ρ [ j..∞] the suffix s j s j+1 . . .. Given a finite word ρ = s 0 s 1 · · · s m , we denote by ρ j the symbol s j for 0 ≤ j ≤ m, and lst(ρ) the symbol s m .
Let AP denote a finite set of atomic propositions. A concurrent game structure (CGS) is a tuple G (S , S 0 , Ag, (Ac i ) i∈Ag , (∼ i ) i∈Ag , (P i ) i∈Ag , ∆, λ), where S is a finite set of states; S 0 ⊆ S is a set of initial states; Ag = {1, ..., n} is a finite set of agents; Ac i is a finite set of local actions of agent i; ∼ i ⊆ S × S is an epistemic accessibility relation (an equivalence relation); P i : S → 2 Ac i is a protocol function such that P i (s) = P i (s ′ ) for every s ∼ i s ′ ; ∆ : S × Ac → S is a transition function with Ac = i∈Ag Ac i being a set of joint actions; and λ : S → 2 AP is a labeling function which assigns each state with a set of atomic propositions. Given a joint action a = a 1 , ..., a n ∈ Ac, we use a(i) to denote the local action of agent i in a.
Intuitively, i (resp. I) signals that agents can only observe partial information characterized via epistemic accessibility relations (resp. complete information with all epistemic accessibility relations being the identity relation), while r (resp. R) signals that agents can make decision based on the current observation (resp. the whole history of observations). We will, by slightly abusing notation, extend Ir-strategies and ir-strategies θ i to the domain S + such that for all ρ ∈ S + , θ i (ρ) = θ i (lst(ρ)). We denote by Θ σ i for σ ∈ {Ir, IR, ir, iR} the set of σ-strategies for agent i.

Outcomes.
A collective σ-strategy of a set of agents A is a function υ σ A assigning each agent i ∈ A with a σ-strategy υ σ A (i) ∈ Θ σ i . For i ∈ A and ρ ∈ S + , we denote the local action υ σ A (i)(ρ) of agent i by υ σ A (i, ρ), and the set Ag \ A by A. Given a state s, two collective σ/σ ′ -strategies υ σ A and υ σ ′ A yield a path ρ, denoted by play(s, υ σ For every state s ∈ S and collective σ-strategy υ σ A of A, the outcome function is defined as follows: i.e., the set of all possible plays that may occur when each agent i ∈ A enforces its σ-strategy υ σ A (i) from the state s. The subscript G is dropped in O σ G when it is clear from the context.

Alternating-Time Temporal Logics
The alternating-time temporal logic ATL/ATL * is an extension of the branching-time logic CTL/CTL * by replacing the existential path quantifiers E with coalition modalities A [Alur et al., 2002]. Intuitively, the formula A φ expresses that the set of agents A has a collective strategy to achieve the goal φ no matter which strategies the agents in A choose. Formally, ATL * is defined by the following grammar: where ϕ (resp. φ) denotes state (resp. path) formulae, q ∈ AP and A ⊆ Ag.
The derived operators are defined as usual: An LTL formula is an ATL * path formula by restricting ϕ to atomic propositions.
The semantics of ATL * is traditionally defined over CGS. When agents's abilities are considered, it is often parameterized with a strategy type σ ∈ T str , denoted by ATL * σ [Bulling and Jamroga, 2014]. Formally, let G be a CGS and s be a state of G, the semantics of ATL * σ (i.e. the satisfaction relation) is defined inductively as follows: (where ̺ is a state s or path ρ) • G, s | = σ q iff q ∈ λ(s); • G, s | = σ A φ iff there exists a collective σ-strategy υ σ A of agents A such that ∀ρ ∈ O σ (s, υ σ A ), G, ρ | = σ φ; • G, ρ | = σ ϕ iff G, ρ 0 | = σ ϕ; Vanilla ATL. ATL is a sublogic of ATL * where each occurrence of the coalition modality A is immediately followed by a temporal operator. Formally, ATL is defined by the following grammar: Remark that the operator R cannot be defined using other operators in ATL with imperfect information [Laroussinie et al., 2008], so is included for completeness.
Given an ATL * formula ϕ, a CGS G and a strategy type σ ∈ T str , the model-checking problem is to determine whether G, s | = σ ϕ or not, for each initial state s of the CGS G.
Some Semantic Issues. We observe that the semantics of ATL/ATL * refers to the agents' abilities in an implicit manner. For the formula A ϕ, the specified σ-strategies only apply to agents in A while the agents in A could still choose beyond σ-strategies (e.g. IR-strategies). In other words, A has a collective σ-strategy to achieve ϕ no matter what the other agents do. When σ is IR as in the original work by [Alur et al., 2002], this interpretation of A ϕ is plausible, as "no matter what the other agents do" is effectively the same as "no matter which strategies the other agents choose". However, when σ is set to be more restricted than IR, agents not in A are still allowed to use IR-strategies.
As mentioned in the introduction, this results in a few shortcomings. From a modeling perspective, arguably agents' abilities should be decided by the practical scenario. Namely, they should be fixed when the model is built, and all agents stick to their respective abilities independent of logic formulae. More concretely, from the semantic perspective, the existing semantics does not take into account the abilities of agents who are not in A, and neglects the (multi-player) game nature in the evolution of MAS. As a result, it may exhibit some counter-intuitive semantics. For instance, consider two formulae A φ and A ′ φ ′ such that agent i ∈ A \ A ′ , i may have different abilities to achieve φ and φ ′ . Let us consider an autonomous road vehicle scenario to see why this is not ideal. There are several autonomous cars which can only observe partial information and have bounded memory. An MAS model G consists of agents in A modeling autonomous cars, and an additional environment agent e. We can reasonably assume that all the car agents use ir-strategies, while e uses IR-strategies. The property A ′ φ expresses that autonomous cars A ′ ⊂ A can cooperatively achieve the goal φ no matter what strategies the other cars and the environment choose. Verifying that G satisfies A ′ φ under the existing semantics would allow car agents A\A ′ to use IR-strategies. If G satisfies A ′ φ, then the result is conclusive, i.e., A ′ φ holds for the system. However, if G invalidates A ′ φ, we cannot deduce that A ′ φ fails because we overestimate the abilities of agents in A \ A ′ when evaluating A ′ φ. In other words, for the formula A ϕ under | = σ where σ IR, it seems to be inappropriate to render the agents not in A extra powers of IR to potentially defeat agents from A when their abilities are actually much weaker.

Agents' Abilities Augmented Concurrent Game Structures
In this section, we introduce agents' abilities augmented concurrent game structures (ACGS in short), which explicitly equip each agent with a strategy type from T str . As such, agents have fixed abilities throughout their lives for a given CGS. Formally, an ACGS is a pair M (G, π), where G is a CGS and π : Ag → T str is a function that assigns a strategy type π(i) to the agent i. We assume that, for each agent i ∈ Ag with π(i) ∈ {IR, Ir}, ∼ i is an identity relation, as agents with perfect information should be able to distinguish two distinct states. Paths of M are defined the same as for the CGS G, but strategies and outcomes of M have to be redefined as follows.

Strategies and Outcomes. Let A be a set of agents. A collective strategy of
Given a state s ∈ S , for every collective strategy ξ A of the set A of agents, the outcome O M (s, ξ A ) of M is the set of all possible paths that may occur when each agent i ∈ A enforces its π(i)-strategy ξ A (i) from state s, and other agents i ∈ A can only choose π(i)-strategies instead of general IR-strategies.
Semantics of ATL/ATL * . The difference of outcomes between ACGS and CGS induces a distinct semantics of ATL * on ACGS than CGS. Let M be an ACGS and s be a state in M, the semantics of ATL * on M is defined similar to the one on CGS, except that the semantics of the state formulae of the form A φ is defined as follows: Remark that this semantics takes into account whether the agents from A have perfect or imperfect information/recall. Given an ACGS M and an ATL/ATL * formula ϕ, the model-checking problem is to determine whether M, s | = ϕ holds, for every initial state s of G. Given a state formula ϕ, let ϕ M denote the set of all the states of M that satisfies ϕ.
The semantics of ATL/ATL * defined on ACGS is different from the one defined on CGS, hence they are incomparable.
Proposition 1 There are an ACGS M = (G, π), an ATL/ATL * formula A φ, and a type σ ∈ T str such that π(i) = σ for all i ∈ A and M, s | = A φ holds, but G, s | = σ A φ.

Model-Checking Algorithms
It has been shown that the Turing Halting problem can be reduced to the model-checking problem of CGS against the ATL formula ϕ = {1, 2} G ok under the iR setting [Dima and Tiplea, 2011], where ok is an atomic proposition. By adapting the proof, we get that: By Theorem 1, we focus on the model-checking problem of ACGS by restricting the function π to Ag → T str \ {iR}. In general, we propose model-checking algorithms which iteratively compute the set of states satisfying state formulae from the innermost subformulae. The main challenge is to compute A φ M . To this end, we first show how to compute A φ M for a simple formula A φ, and then present the more general algorithm. An ATL/ATL * formula A φ is called simple if φ is an LTL formula.
Let us fix an ACGS M = (G, π) with G = (S , S 0 , Ag, (Ac i ) i∈Ag , (∼ i ) i∈Ag , (P i ) i∈Ag , ∆, λ) and a simple formula A φ. Given a set of agents A ′ and a strategy type σ ∈ T str , we denote by A ′ σ the set {i ∈ A ′ | π(i) = σ}.

Model-Checking Simple ATL Formulae
For a simple ATL formula A φ, it is easy to see that whether agents in A have perfect call or not does not matter if these agents have perfect information abilities.
Proposition 2 Given an ACGS (G, π) with π : Ag → T str \ {iR}, and a simple ATL formula A ϕ, let π ′ be a function such that for By Proposition 2, all the agents in A with IR-strategies can be seen as agents with Ir-strategies. All the agents in Ag with Ir-strategies can be seen as agents with ir-strategies (i.e., all the epistemic accessibility relations of them are the identity relation). Therefore, we can assume that π(i) = ir for all i ∈ A, and π(i) ∈ {ir, IR} for all i ∈ A. For two collective strategies ξ A and ξ A ir , let M(ξ A , ξ A ir ) = (G ′ , π) be the ACGS obtained from (G, π) by enforcing strategies ξ A and ξ A ir , namely, by removing transitions whose actions of agents in A ∪ A ir do not conform to ξ A and ξ A ir . We have that Computing ∅ φ M(ξ A ,ξ A ir ) amounts to CTL modelchecking, which can be done in polynomial time (and thus et al., 1983]. Since the number of strategies ξ A and ξ A ir is finite, we get that:

Model-Checking Simple ATL * Formulae
We compute A φ M by a reduction to the problem of computing the winning region of a turn-based two-player parity game. We first introduce some basic concepts which will be used in our reduction.
Let inf(ρ) be the set of states visited infinitely often in ρ. An infinite word is recognized by A if A has a run ρ over this word such that min p∈inf(ρ) R(p) is even. For the LTL formula φ, one can construct a DPA A φ = (P, 2 AP , δ, p 0 , R) with 2 2 O(|φ|) states and rank k = 2 O(|φ|) such that A φ recognizes all the ω-words satisfying φ [Piterman, 2006].
A (turned-based, two-player) parity game P is a tuple ( Given a strategy θ 0 for Player-0 and a strategy θ 1 for Player-1, let G θ 0 ,θ 1 be the play where Player-0 and Player-1 enforce their strategies θ 0 and θ 1 . θ 0 is a winning strategy for Player-0 if min s∈inf(G θ 0 ,θ 1 ) Ξ(s) is even for every strategy θ 1 of Player-1. The winning region of Player-0, denoted by WR 0 , is the set of vertices from which Player-0 has a winning strategy.
A partial strategy of A ′ is a partial function f : and g ⊥ ∈ G ir denote the partial strategy such that dom(g ⊥ ) = ∅. Given a state s, let F s IR be the set of functions f : A IR → i∈A IR Ac i such that for every i ∈ A IR , f (i) ∈ P i (s). Let F IR := s∈S F s IR and Π ir : We define a parity game k} is a function such that for every s ∈ S , ir , ∀ f ∈ F IR and ∀G ∈ Π ir , E is defined as follows: where G ′ ⊆ Π ir is the largest set such that the follows hold: for every g ′ ∈ G ′ , In this reduction, f ⊤ encodes a collective ir-strategy of agents in A ir , the collection of f 's in one play of P φ encodes a collective IR-strategy of agents in A IR , and each g ∈ G encodes a collective ir-strategy of agents in A ir . The imperfect information abilities of agents are ensured by the definitions partial strategies.
To check whether s ∈ A φ M , P φ starts with the vertex s. At first step, Player-0 chooses a function f ⊤ ∈ F ⊤ ir meaning that all agents in A ir choose an ir-strategy. Next, P φ moves from (s, f ⊤ ) to (s, p 0 , f ⊤ , {g ⊥ }) which let the DPA A φ start with p 0 (note that Player-1 has only one choice at this step). At a vertex (s, p, f ⊤ , G) controlled by Player-0, Player-0 chooses actions for agents in A IR by choosing one function f ∈ F s IR . Then Player-1 chooses actions for agents in A with respect to the chosen actions of agents in A ir tracked by G. These selections of actions together with f ⊤ and G determine a joint action a, based on which P φ moves to (s ′ , δ(p, λ(s)), f ⊤ , G ′ ) such that s ′ is the next state of s under a, δ(p, λ(s)) is the next state of p in A φ which allows to mimics the run of A φ over the ω-word induced by the play of M. During this step, f is dropped from the vertex of P φ , as f corresponds to actions of agents in A IR and needs not to track. The actions of agents in A ir are tracked by computing G ′ from G. This ensures imperfect recall abilities of agents in A ir . We then can get that: The winning region of Player-0 in P φ can be computed in polynomial time of |V| · |E| · 2 k [Jurdzinski, 2000]. In this reduction, each G contributes at most O(|S |) sets of G ′ . Therefore, |V| · |E| is exponential in |G| · 2 |φ| . Recall that k = 2 O(|φ|) .

Lemma 3 For the simple ATL * formula A φ,
A φ M can be computed in 2EXPTIME.

12
Compute A φ M by Lemma 1 or 3; 13 return A φ M ;

The Overall Algorithm
We now present the overall procedure, which computes ϕ M from the innermost subformulae. Algorithm 1 shows the pseudo code, which takes an ACGS M = (G, π) and an ATL/ATL * formula ϕ as inputs, and outputs ϕ M . Then M satisfies ϕ iff the set of initial states of M is a subset of ϕ M . We also incorporate epistemic modalities Hoek and Wooldridge, 2003;Cermák et al., 2014] into our algorithm with the following semantics: . K i ϕ, E A ϕ, D A ϕ and C A ϕ respectively denote that i knows, every agent in A knows, agents in A have common knowledge and agents in A have distributed knowledge, on the fact ϕ. The ATL (resp. ATL * ) logic extended with these epistemic modalities is called ATLK (resp. ATLK * ) logic. Given a state s ∈ S and a binary relation ⋍⊆ S × S , we denote by [s] ⋍ the set {s ′ ∈ S | s ⋍ s ′ }.
By Lemma 1 and Lemma 3, the model-checking problem for ATLK and ATLK * on ACGS can be solved in PSPACE and 2EXPTIME respectively. As the model-checking problem of ATL * on CGS under IR setting is 2EXPTIMEcomplete [Alur et al., 2002], we have that Theorem 2 The model-checking problem for ATLK * (resp. ATLK) on ACGS is 2EXPTIME-complete (resp. in PSPACE).  In this CG game, there are several agents modelling workers and an environment agent. Each worker works for the benefit of a castle, and the environment keeps track of the Health Points (HP) of castles. Each castle preserves an HP valued up to 3, and 0 means it's defeated. Workers are able to attack a castle which they don't work for, or defend the castle which they work for, or do nothing. The castle gets damaged if the number of attackers is greater than the number of defenders, and the differences influence its HP. In this model, the number of states is 8000 × 4n, the environment agent has 1 local action, and each worker agent has 4 local actions, where n denotes the number of workers.

Implementation and Experiments
In this experiment, we consider ACGS consisting of three worker agents w 1 , w 2 , w 3 and an environment agent e, where worker w i works for the castle c i .
• ϕ 1 ≡ {w 1 , w 2 } F(castle3De f eated): expresses that workers w 1 and w 2 can make castle c 3 defeated, no matter which strategies other agents use.
• ϕ 2 ≡ {w 1 , w 2 } F(allDe f eated): expresses that workers w 1 and w 2 can make all the castles defeated, no matter which strategies other agents use.
We observe that: (1) the strategy types of agents do affect the performance and results. In particular, the time significantly increases when w 3 is ir-typed while w 1 or w 2 is IRtyped; (2) Lemma 1 is more efficient when both w 1 and w 2 are ir-typed; otherwise Lemma 3 is more efficient. This is because the number of possible strategies of w 1 and w 2 is small (using Lemma 1) if both w 1 and w 2 are ir-typed.

Related Work
The family of alternating-time temporal logics ATL, ATL * and alternating µ-calculus (AMC) [Alur et al., 2002] for rea-soning about games was introduced with motivations partially from MAS. Model-checking algorithms were also given with IR-strategies. [van der Hoek and Wooldridge, 2003] extended ATL with knowledge operators and proposed corresponding model-checking algorithms. In their work, epistemic accessibility relations are considered in the interpretation of knowledge operators, but not for the strategies and outcomes. This means that agents still use IR strategies for coalition modalities A ψ. This issue was discussed in [Jamroga, 2003] which proposed an idea of iRstrategies. [Schobbens, 2004] introduced the notion of imperfect recall into ATL/ATL * , and investigated the modelchecking problem for ATL/ATL * under four different strategic types. Importantly, with iR-strategies the model-checking problem becomes undecidable [Dima and Tiplea, 2011].
[Bulling and Jamroga, 2011] studied the semantics of AMC and proposed a model-checking algorithm for the alternationfree fragment under the imperfect information setting.
[Bulling and Jamroga, 2014] further conducted a comprehensive comparison of variants of ATL/ATL * with different strategic types. The study corroborate that the agents' abilities play a prominent role in logic semantics.
In the previous work, strategies are revocable, i.e., when it comes to achieve a goal in the (nested) sub-formula, previously selected strategies are deleted. [Ågotnes et al., 2007] introduced a variant of ATL with irrevocable strategies under the imperfect recall setting. It was generalized into ATL/ATL * with strategy contexts [Laroussinie and Markey, 2015], which allows agents to drop or inherit previously selected strategies.
Two strategic logics were introduced by [Chatterjee et al., 2010] and [Mogavero et al., 2014] and the model-checking problem was investigated therein under the IR-setting. Strategic logics extend LTL with first-order quantifications over strategies which naturally captures the multi-player game nature in the evolution of MAS. [Cermák et al., 2014] introduced knowledge operators in the strategic logic of [Mogavero et al., 2014] and proposed a model-checking algorithm with ir-strategies. Here all agents must take ir-strategies (so the potential inconsistency can be ruled out), but no other abilities are considered. To gain decidability under iR-setting, specific restrictions on the abilities of the agents were proposed [Berthon et al., 2017a;Berthon et al., 2017b;Belardinelli et al., 2017a;Belardinelli et al., 2018].
Our work is orthogonal to the existing work which defines the agents' abilities at the semantics level, but takes a more syntactic level by strengthening the model.

Conclusion and Future Work
In this paper, we discussed the problem of existing semantics of ATL/ATL * , and advocated the approach to make agents' abilities explicit in modeling. For this purpose, we introduced an extension of standard CGS model, i.e. ACGS, which define agents' abilities at the syntactic level of the system model. We explored the effects of strategy types in the semantics, in particular model-checking, of ATL/ ATL * over ACGS, and pro-vided model-checking algorithms with identified complexity. The algorithms are implemented in a tool MCMAS-ACGS, which has been applied to several applications to demonstrate the feasibility of the approach. This work represents the first systematic study towards different agents' abilities at the syntactic level, which is in contrast to the previous approaches at the semantic level.
Currently we use ATL/ATL * as the specification, but the methodology can be extended to other logics such as Strategy Logic, and other agents' abilities such as strategy contexts. Several questions are left open such as axiomatization and satisfiability problem. We leave them for future work.

A Proof of Theorem 1
Given a Turing machine M, let G = (S , {s 0 }, Ag, (Ac i ) i∈Ag , (∼ i ) i∈Ag , (P i ) i∈Ag , ∆, λ) be the CGS constructed as in [Dima and Tiplea, 2011] and s 0 be an initial state in G such that G, s 0 | = iR {1, 2} G ok iff M does not halt on the empty word. Let M = (G, π) be an ACGS such that for every agent i in M, π(i) = iR if i ∈ {1, 2}, otherwise π(i) = IR. Therefore, The result immediately follows.

B Proof of Proposition 2
Without loss of generality, we assume that for every i ∈ Ag, π(i) Ir, as Ir can be seen as a special case of ir. We first construct the tree-unfolding of M from the state s.
(⇐) Suppose M * s , s | = A ϕ, let ξ * A be the collective strategy such that for every path ρ ∈ O M * s (s, ξ * A ), M * s , ρ | = ϕ. We assume that there is a total order on set S + , and denote by min(U) the minimal one of the set of states U ⊆ S + with respect to the order . Let ξ A be the function such that for every i ∈ A and s ′ ∈ S , ξ A (i)(s ′ ) = ξ * A (i)(min({ρ ∈ S + | lst(ρ) = s ′ })). First, we show that ξ A is a collective strategy of A in M ′ . Consider an agent i ∈ A and two states s 1 , s 2 ∈ S such that s 1 ∼ i s 2 , if π(i) IR, then for each pair of states ρ 1 , ρ 2 ∈ S + such that lst(ρ 1 ) = s 1 and lst(ρ 2 ) = s 2 , we have ρ 1 ∼ * i ρ 2 , which implies that ξ * . The result immediately follows from the fact that λ * (ρ) = λ(lst(ρ)) for every ρ ∈ S + .
⊓ ⊔ Recalling that [Alur et al., 2002] observed that both semantics of ATL under Ir-strategies and IR-strategies coincide for CGS . This result was generalized and formally proved in infinite CGS (i.e., no finiteness with respect to the set of states and actions) (cf. Proposition 1 [Bulling and Jamroga, 2014]). Proposition 2 can be seen as a generalization of the result of [Alur et al., 2002] and could be extended to the infinite ACGS similar to [Bulling and Jamroga, 2014].

C Effects of Strategy Types
By restricting all the strategy types to IR, straightforwardly we have: Proposition 3 Let M = (G, π) be an ACGS where for each i ∈ Ag, π(i) = IR. For each state s of M and ATL * formula ϕ, G, s | = IR ϕ iff M, s | = ϕ.
Proof. By applying structural induction, it suffices to show that the result holds for formulae of the form A φ. By applying induction hypothesis, for every path ρ, the following holds: . Each i ∈ A has same sets of possible IR-strategies in G and M, hence ⊓ ⊔ In light of Proposition 1 and Proposition 3, we shall investigate the effects of strategy types by considering ACGS with various different setups of strategy types.
Lemma 4 Let A be a set of agents and s be a state of a CGS G. For two functions π 1 , π 2 : Ag → T str with π 1 A π 2 , and any collective strategy ξ A of A, we have: Lemma 4 reveals the effect of strategy types of A on the outcomes. It is easy to observe that if π 2 (i) = σ for all i ∈ A, then for every collective σ-strategy υ σ (3) negations ¬ only appear in front of atomic propositions. For example, A X q is positive, while ¬ A X q is not positive. Given a formula ϕ, let Ag ϕ denote the set of agents that appear in ϕ. By Lemma 4, we have: Proposition 4 Let G be a CGS, s be a state of G and ϕ be a positive ATL/ATL * formula. For π 1 , π 2 : Ag → T str such that π 1 Ag ϕ π 2 , if (G, π 2 ), s | = ϕ, then (G, π 1 ), s | = ϕ.
Proof. By applying structural induction, it suffices to show that the result holds for formulae of the form A φ. We suppose (G, π 2 ), s | = A φ, otherwise Item 1 immediately holds.
⊓ ⊔ For positive ATL/ATL * formulae ϕ, even if the agents of Ag ϕ have the same strategy types in ACGS (G, π) and CGS G, verifying G against ϕ under σ will examine more behavior than verifying (G, π) against ϕ, where i ∈ Ag ϕ and π(i) = σ. Therefore, if the behavior of a MAS is exactly modeled as an ACGS M rather than a CGS G with strategy type σ, verifying G against ϕ under σ may lead to incorrect result. However, more restrictions on strategy types and ATL/ATL * formulae can make them coincide, as the following proposition shows.

D.1 Experiment on Dining Cryptographers Protocol
Dining cryptographers protocol is one of anonymity protocols aimed at establishing the privacy of principals during an exchange . The dining cryptographers protocol can be modeled as MAS. In this game, n cryptographers share a meal around a circular table. Either one of them or their employer paid for the meal. They are curious whether it was sponsored by their employer without revealing the identity of the payer (if one of them did pay). The protocol works as follows: each cryptographer 1) tosses a coin and shows the outcome to his right-hand neighbour, 2) announces whether the two coins agree or not if he/she is not payer, otherwise announces the opposite of what he/she sees. Their employer is the payer if an even number of cryptographers claiming that the two coins are different, otherwise not. For experimental purpose, we allow the cryptographer who paid for the meal announces either the two coins agree or not no matter what he/she saw.
In this experiment, n ranges from 3 to 10, a cryptographer whose is payer and one of other cryptographers use ir-strategies, the others all use IR-strategies. We verify three formulae ψ 1 , ψ 2 and ψ 3 , where ψ i expresses that if the number of "saydifferent"  is odd and the i-th cryptographer is not the payer, then he/she knows that the bill is paid by one of the others, but cannot tell exactly who is the payer. For instance, in the three cryptographers case, The results are shown in Table 2, where column 1 gives the number of cryptographers, column 2 gives the number of states, columns 3-5 (resp. columns 6-7) show the total time of respectively verifying ψ 1 , ψ 2 and ψ 3 using Lemma 1 (resp. Lemma 3). Both ψ 1 and ψ 2 are satisfied by all the models, while ψ 3 not. We observe that Lemma 1 is more efficient than Lemma 3, as the coalitions in all the formulae are ∅. From this experiment, one may conclude the reasonable scalability of our tool.

D.2 Experiment on Book Store Scenario (BSS)
The BSS model depicts a deal between two agents: a supplier (S) and a purchaser (P) . Initially, S is waiting for an order from P, and P is ready for initiating a trade. Upon receiving an order of some e-good from P, S can make a decision whether accepts the order or not do, and later notify P. If S accepts, then P can pay fee. Once paid, S can either reject the payment or accept and deliver the good. If P received the good, then trade is completed. During the trade, P can revoke the order, both S and P can terminate the trade, after which the information of the trade should be symmetric at any time. In this model, the S has 15 local states and 13 actions, and P has 12 local states and 7 actions. In this experiment, we verify the model against the following formulae. ϕ 1 ≡ ∅ G((S &P no T ) → (K S {S , P} F trade end)) expresses that "if neither S nor P terminates the trade (i.e., S &P no T is ture), then S knows that they can cooperatively complete the trade eventually" always holds. ϕ 2 ≡ {S , P} (S &P no T U(trade end ∧ ¬trade success)) expresses that the trade can ends by P asking for refund.
The results are shown in Table 3. Each row presents the result of one of strategy type combinations of S and P, for instance, (IR, ir) denotes that S has an IR-strategy while P has ir-strategy. Columns 2-4 (resp. columns 5-7) show total time and result of verifying ϕ 1 (resp. ϕ 2 ) using Lemma 1 (resp. Lemma 3). The results of ϕ 1 confirm that strategy types affect the truth of formulae. Lemma 1 performs better than Lemma 3 both on ϕ 1 and ϕ 2 in this experiment.


extfieldsofstudy:
Physics
Computer Science

text:
Using “GUI” programming environments to calculate the level of service at intersections

Inside busy cities, road congestion is common due to the large number of vehicles. Modifying existing infrastructure only for certain periods of peak time or for specific events involves high costs, so it is desirable to use its full capacity. At the same time those in charge of traffic optimization level infrastructure and coordination are municipalities. Based on the above premises it can be seen that for a better optimization of traffic in a critical situation and require the use of unique software for calculating traffic light times, service level, etc. The paper aims at using software platforms in “GUI” applications to be used in the context of dynamic changes in road traffic.

Introduction
The level of service indicates the degree of use (load intersection, through the transit of vehicles and a degree of creating congestion) at an intersection, translated in time constant [1] [2].
These criteria of analysis are necessary for rapid examinations imposed a junction and offers some quick values that can identify in real time where the intersection and the problems studied.
The maximum service flow rate is defined as the maximum number of vehicles that can cross a uniform road section over a given time period (preferably 15 minutes) under certain road, traffic and control conditions to maintain a certain level of service projected [8].
They are highlighted as follows: • Level A: characterizes free flow, with individual users, virtually unaffected by the presence of other traffic vehicles. • Level B: characterizes stable flows with a high degree of freedom in choosing speed and operating conditions that have little influence on other road users. • Level C: characterizes restricted flows that remain stable but interact with other participants in traffic flow. The overall comfort and safety level drops considerably. • Level D: characterizes the high-density flow in which speed and maneuverability are severely restricted, and comfort and safety are low even if the flow remains stable. • Level E: characterizes unstable flow at or near the capacity limitation, with the lowest level of comfort and safety. • Level F: corresponds to conditional traffic in which the number of vehicles arriving at a point exceeds the service possibilities and therefore creates the conditions of tail formation (traffic jams), low comfort and increased accident risk.

Designing the graphical user interface (GUI)
Designing the graphical user interface was done using the MATLAB platform, GUI mode. Figure 2 shows how the buttons, windows, etc. were arranged in the design interface of the user interface, and also figure 2 represents the graphical interface in the running mode, i.e. as seen by the user [3] [4].

Figure 2. Working window for the GUI
The main window ( Figure 2) is divided into four distinct zones [5] in order to enter the input data for the level of service calculation:  Area 1 is the area where inputs will be entered for vehicle volumes for each lane of bands existing on that arm.  Zone 2 represents the area where the saturation flow for each group of bands will be entered.  Zone 3 is the area where all the values associated with the types of traffic lights are to be entered, depending on the signal plan configuration.  Zone 4 is the group area of the control delay calculation buttons for the band groups, arms and intersection studied. In order for the software not to be stuck in an infinite loop or loop due to errors that may occur when data is entered, codes are provided on all critical areas to provide error windows. For example, one of these errors is found at the introduction of traffic phases, so this window appears if less than two phases are entered [6].

Conceiving the software instruction code
This subchapter presents the main sequence of instructions, sequences that calculate, check, and display the service level values at the intersection studied. Figure 5 shows the code for initializing the graphical interface, entering the text boxes with the default value of 0, and resetting the constants and the database to calculate the service level.  The next step is to retrieves user input into the text boxes and converts them into input variables to be entered into the service level code. After the values are transformed from string to integer variables the code for calculating the traffic capacity that takes into account the intersection geometry and all the traffic volumes stored in the database when the software is opened, is called, as can be seen in the figure 7.
After the traffic capacity is calculated the code for the flow-capacity calculation is called. All variables used in the program code are retrieved from the database, created temporarily by the software. The next code called is for the calculation of the uniform delays and incremental delays at the intersection studied, and uses "for" and "if" instructions to be able to calculate the delays values as accurately as possible depending on the intersection geometry. The next step is the call of the code for the calculation of control delays for every group of movement, the aggregation of the control delay for

Case study for a junction in Craiova city
Verification of the LOS and at the same time the functionality of the computing application made in the MatLab platform will be accomplished with data taken from an intersection in Craiova. The intersection chosen is a transfer area between two neighborhoods and the city center, as well as a critical point thanks to the municipal hospital in Craiova, as can be seen from Figure 8.  In order to be able to check the service level of the intersection studied, it is necessary to provide the software with the data on the intersection traffic lights. Thus, the greening times at the intersection reported for each phase (the northern phase has 32 seconds, the south phase has 29 seconds, the phase is 16 seconds, and the west phase has 15 seconds), were measured also. It can also be seen that the duration of the traffic light cycle is 112 seconds, and the intersection protection and release times at the end of each phase are 3 seconds of yellow and 2 seconds of red. After running the input data through the LOS calculation software, we obtained the specific intersection results studied divided into delays related to movement groups, arm and intersection.

Conclusions
The utility behind the LOS MATLAB application is being designed to help in the following aspects:  A rapid verification tool for intersection delays;  Graphical user interface as simple and intuitive as possible;  Creating the programming procedures so as to use as few resources as possible;  It is an application designed for rapid calculation so that it can be easily used by personnel without special training in the field;  Can be used on the field or in the office by changing input data very quickly to see the differences in delay time between different traffic flows;  The application can rapidly evaluate different types of signal plans based on the LOS results.