Pietro Barbiero

Giovanni Squillero
email: squillero@polito.it

Alberto Tonda
email: alberto.tonda@inrae.fr

Predictable Features Elimination: An Unsupervised Approach to Feature Selection

We propose an unsupervised, model-agnostic, wrapper method for feature selection. We assume that if a feature can be predicted using the others, it adds little information to the problem, and therefore could be removed without impairing the performance of whatever model will be eventually built. The proposed method iteratively identifies and removes predictable, or nearly-predictable, redundant features, allowing to trade-off complexity with expected quality. The approach do not rely on target labels nor values, and the model used to identify predictable features is not related to the final use of the feature set. Therefore, it can be used for supervised, unsupervised, or semi-supervised problems, or even as a safe, pre-processing step to improve the quality of the results of other feature selection techniques. Experimental results against state-of-the-art feature-selection algorithms show satisfying performance on several non-trivial benchmarks.

Introduction

The aim of most Machine Learning (ML) algorithms is to build a predictive model starting from the feature values of a given training set. State-of-the-art algorithms are usually quite effective at tacking problems with huge number of samples, yet they might face issues if the number of features is huge. An increase in the dimensionality of the problem, in fact, may correspond to a much steeper increase of the search space, impairing the optimization of the models, or creating other problems usually described with the vague expression "curse of dimensionality" [START_REF] Barbiero | Modeling generalization in machine learning: A methodological and computational study[END_REF].

A reduction in the number of variables may be obtained either by means of feature extraction or feature selection. Techniques in the former group, such as principal component analysis (PCA) or autoencoders, built a new, more compact set of features out of the original one. Feature selection techniques, on the other hand, aim at finding a subset of the original features, that still allows ML algorithms to build reliable predictive models. There might be practical reasons to opt for the second class of techniques, for instance, if each single feature has a cost because it needs to be physically measured. Moreover, predictive models that use tens or hundreds of features are de facto black boxes, quite hard if not completely impossible to interpret. Identifying the key features involved in a problem can make the final model more human-readable: problems involving genomic data [START_REF] Bermingham | Application of high-dimensional feature selection: Evaluation for genomic prediction in man[END_REF] would greatly benefit from the possibility to find understandable correlations; improvement could be noticeable in the field of visualization [START_REF] Tsai | Dimensionality reduction for computer facial animation[END_REF]; and the general need for explainable artificial intelligence (XAI) is only but increasing in present days.

Starting from the observation that the value of a redundant feature, by definition of redundancy, can be inferred using the information contained in the other features, we propose Predictable Features Elimination (PFE). PFE is an unsupervised, model-agnostic, wrapper approach for feature selection: In a first step, each feature is scored and ranked using a statistical measure; then, starting from the lowest-ranked feature, an auxiliary ML model is trained to predict that feature using all the others; if the performance of the model exceeds a given quality (for example, R2 > 0.95), the information it provides is assumed to be redundant and the feature is removed. The procedure then iterates to the next feature, and once all features have been analyzed in this way, those remaining represent the final subset.

Experimental results on several non-trivial benchmarks from the OpenML repository [START_REF] Vanschoren | Openml: Networked science in machine learning[END_REF] show that the proposed approach is competitive with the state of the art in the field, obtaining feature subsets that are either more informative or smaller than competing feature selection algorithms. The main contributions of this paper are:

-We describe Predictable Features Elimination, a new unsupervised, modelagnostic, wrapper approach for feature selection. -We compare the performance of PFE against state-of-the-art feature selection algorithms, showing the advantages of using our method especially on large datasets. -We compare PFE against other algorithms on an artificial dataset, specifically designed for assessing feature selection algorithms. Results show that the features selected by PFE include almost all meaningful information.

Feature Selection in Machine Learning

Feature Selection (FS) is the process of identifying the features of a data set in order to obtain a minimal, informative subset. Features may not be part of this subset for two main reasons: they might be unrelated to the underlying nature of the problem, just adding noise; they might be heavily correlated with others features, adding no relevant information for the task. Applications range from face recognition [START_REF] Vignolo | Feature selection for face recognition based on multi-objective evolutionary wrappers[END_REF] to medicine [START_REF] Zhou | Multi-objective based radiomic feature selection for lesion malignancy classification[END_REF], while approaches can be divided into two categories [START_REF] Guyon | An introduction to variable and feature selection[END_REF]: filters that score features according to a criterion (often a statistical test); and recursive procedures (forward or backwards) that attempt to reduce the features to a small set of non-redundant ones [START_REF] Lewis | The characteristic selection problem in recognition systems[END_REF][START_REF] Chien | On the generalized karhunen-loève expansion (corresp.)[END_REF].

Several different metrics that have been proposed to assess the content of information of a given feature subset: from mutual information [START_REF] Kozachenko | Sample estimate of the entropy of a random vector[END_REF] to analysis of variance [START_REF] Fisher | Xv.-the correlation between relatives on the supposition of mendelian inheritance[END_REF]. However, only some of the information content of a feature subset can be assessed through such metrics, as taking into account the contribution of non-linear combinations of features would be too computationally expensive. Recursive procedures on the other hand frequently rely on more complex measurements, usually a goodness-of-fit [START_REF] Guyon | An introduction to variable and feature selection[END_REF] for a model wrapped inside the feature selection, sometimes combined with regularization terms [START_REF] Weston | Feature selection for svms[END_REF]. In some cases, the number of features and the quality of the fit are evaluated separately, and each candidate subset is placed on a Pareto front [START_REF] Xue | A survey on evolutionary computation approaches to feature selection[END_REF].

Literature reports FS techniques as complex as exploiting evolutionary algorithms (EAs) [START_REF] Cilia | Variable-length representation for ec-based feature selection in high-dimensional data[END_REF][START_REF] Xue | A survey on evolutionary computation approaches to feature selection[END_REF], with single-or multi-objective approaches [START_REF] Hamdani | Multi-objective feature selection with nsga ii[END_REF][START_REF] Xue | Multi-objective feature selection in classification: a differential evolution approach[END_REF][START_REF] Barbiero | A novel outlook on feature selection as a multi-objective problem[END_REF]. Anyhow, the most popular approach in literature is still probably the 20-year old Recursive Feature Elimination (RFE) [START_REF] Guyon | Gene selection for cancer classification using support vector machines[END_REF], a supervised, wrapper methodology that iteratively removes the worst features based on the performance of the target model.

Predictable feature elimination

The algorithm Predictable Features Elimination, the approach we present in this work, stems from the observation that if the distribution of a feature f r can be approximated by using the information of other features, then f r is likely to be almost redundant and of little importance for whatever model will be eventually built. Algorithm 1 summarizes the main steps of the training process.

PFE requires the user to provide two parameters: an auxiliary machine learning model g and a threshold σ. The first is the model used to discriminate non-redundant features, the second, the acceptable loss of information. A suboptimal choice of g would cause some features to be erroneously marked as non-redundant, increasing the size of the feature set, but probably not affecting the quality of the final model. On the other hand, a low σ is likely to make PFE select a very small set of features, but also to impair the quality of the final model.

The algorithms is composed of two phases: an initialization and the main loop. In the former, features are ranked according to their mutual average linear correlation. First, the feature correlation matrix C is computed:

C = diag(K XX ) 1/2 K XX diag(K XX ) 1/2 (1)
where K XX is the auto-covariance matrix of the input matrix X:

K XX = E[(X -E[X])(X -E[X]) T )] (2)
The correlation matrix C is used to estimate the amount of mutual information shared among features. By summing up the rows of C, we obtain for each feature an approximation of the amount of information which can be obtained using all the other features:

κ = i C(i, •) (3)
The more a feature is correlated with others, the lower the chances that the feature may contain exclusively useful information. Hence, by ranking features according to their mutual average linear correlation, we will obtain an ordered list of their significance:

κ s = sort(κ) (4)
Starting from the feature with the highest rank f , the ML model g is trained on the remaining features using the f -th feature as a target variable y:

y = X(•, κ s (f )) (5)
The performance of g is assessed on a validation set X val using the coefficient of determination R2. If the validation score is greater than the user defined threshold σ, then it means that the model g represents an accurate nonlinear association between the f -th feature and the other features. Hence, the chances that the f -th feature may contain exclusively useful information are low. Therefore, the f -th feature should be safely removed from the feature set and the process may continue using the following feature in the ranking. The algorithm stops when more than half of the features have been analyzed.

Algorithm 1: Predictable feature elimination

Input: data X ∈ R n,d , model g, threshold σ ∈ [0, 1] Initialize C = corr(X T ) Initialize κ = i C(i, •) Initialize κs = sort(κ) for f = 1 to d/2 do Initialize y = X(•,

Theoretical foundation

The following theorems yield a theoretical justification for the proposed approach. Besides, they show how the performance loss can be formally estimated by providing upper bounds in worst case scenarios.

Theorem 1 (Elimination for linear combinations). Let F be a set of features F = {f 1 , . . . , f d } and F be a subset of F such that f i / ∈ F . If the feature f i is a linear combination of the feature set F , then fitting a linear model using F is equivalent to fitting the same model using F ∪ {f i }.

Proof. By definition f i is a linear combination of F , hence:

f i = j =i w j f j (6)
which can be written as:

f i = F w (7)
where F ∈ R n×d is a matrix whose columns are features in F and w ∈ R d is a row vector containing the weights of the linear combination.

A linear model g trained using the matrix F can be written as:

y = g(F ) = F w g = j∈d f j w g j (8)
Let F be the matrix whose columns correspond to the features in F ∪ {f i }, then the model g trained on F can be written as:

y = g(F ) = F w g = k∈d f k w g k = j∈d f j w g j + w g i f i = j∈d f j w g j + w g i j∈d f j w j = j∈d f j w g j + j∈d f j w j w g i = j∈d f j w g j w j w g i = j∈d f j ω g j ( 9
)
Theorem 2 (Approximate elimination for linear combinations). Let F be a set of features F = {f 1 , . . . , f d } and F be a subset of F such that f i / ∈ F . If the feature f i can be written as a linear combination of the feature set F with an additional term , then the upper bound of the training error obtained by fitting a linear model on F instead of F ∪ {f i } is at most . Proof. By definition f i can be written as a linear combination of the feature set F with an additional term , hence:

f i = j =i w j f j + ( 10
)
which can be written as:

f i = F w + ( 11
)
where F ∈ R n×d is a matrix whose columns are features in F and w ∈ R d is a row vector containing the weights of the linear combination.

A linear model g trained using the matrix F can be written as:

y = g(F ) = F w g = j∈d f j w g j ( 12
)
Let F be the matrix whose columns correspond to the features in F ∪ {f i }, then the model g trained on F can be written as:

y = g(F ) = F w g = k∈d f k w g k = j∈d f j w g j + w g i f i = j∈d f j w g j + w g i j∈d f j w j + = j∈d f j w g j + j∈d f j w j w g i + = j∈d f j w g j w j w g i + = j∈d f j ω g j + (13)
Theorem 3 (Approximate elimination). Let F be a set of features F = {f 1 , . . . , f d } and F be a subset of F such that f i / ∈ F . Let g and h be two nonlinear models with equivalent capacity. If the feature f i can be written as a function of F through h with an error term , then the upper bound of the training error obtained by fitting g on F instead of F ∪ {f i } is at most η(g, ).

Proof. By definition f i can be written as a nonlinear function h of the feature set F with an additional term , hence:

f i = h(F ) + ( 14
)
where F ∈ R n×d is a matrix whose columns are features in F . A nonlinear model g trained using the matrix F can be written as:

y = g(F ) (15)
Let F be the matrix whose columns correspond to the features in F ∪ {f i }, then the model g trained on F can be written as:

y = g(F ) = g(F , f i ) = g(F , h(F ) + ) (16)
Since the information in h(F ) can be obtained from F by applying the function h, then g can be fitted on F without information loss by discarding h(F ):

y = g(F ) = g(F , ) = g(F ) + η(g, ) (17)
where η is a function of g and .

The following definition can be used to monitor in an unsupervised way the performance loss when multiple features are recursively eliminated. In some applications, this lemma may be used to derive alternative stopping conditions. Definition 1 (Validity of feature elimination). Let λ be an upper bound of the performance loss required for a specific application and let {η 1 , . . . , η k } be a sequence of training errors obtained by performing k steps of feature elimination. The sequence of k feature elimination steps is valid if and only if λ ≤ k i=1 η i .

Experimental Evaluation

The proposed approach has been implemented from scratch in Python 3, using only open-source libraries [START_REF] Pedregosa | Scikit-learn: Machine learning in python[END_REF][START_REF] Li | Feature selection: A data perspective[END_REF]; the source code, including all the parameter values used in the experiments, is available under the European Union Public Licence (EUPL) from GitHub 4 .

All experiments have been run on the same machine equipped with an AMD EPYC 7301 16-core processor running at 2 GHz, and with 64 GiB memory.

Experimental setup

The performance of predictable feature elimination is compared over a 10-fold cross-validation against both supervised and unsupervised feature selection algorithms [START_REF] Pedregosa | Scikit-learn: Machine learning in python[END_REF][START_REF] Li | Feature selection: A data perspective[END_REF]: Laplacian score for feature selection (lap score [START_REF] He | Laplacian score for feature selection[END_REF]), spectral feature selection for unsupervised clustering (SPEC [START_REF] Zhao | Spectral feature selection for supervised and unsupervised learning[END_REF]), multi-cluster feature selection (MCFS [START_REF] Cai | Unsupervised feature selection for multi-cluster data[END_REF]), non-negative spectral feature selection (NDFS [START_REF] Li | Unsupervised feature selection using nonnegative spectral analysis[END_REF]), regularised discriminative feature selection (UDFS [START_REF] Yang | L2, 1-norm regularized discriminative feature selection for unsupervised[END_REF]), and recursive feature elimination (RFE [START_REF] Guyon | Gene selection for cancer classification using support vector machines[END_REF]).

A Ridge classifier has been used both as the internal estimator for RFE, and to discard redundant features in PFE. It must be noted that PFE is agnostic to the choice of the estimator g, as far it is not significantly superior to the one eventually used in the final model -the underlying assumption being that if a feature can be predicted by g, it may as well be inferred by the final model.

Ridge has been chosen both for its high train speed and its generalization ability in a variety of experimental settings.

For all the experiments the σ parameter has been set to the default 0.9, a reasonable value that leads to the removal of at most half of the original features. For the sake of comparison, all the other feature selection algorithms are set in order to provide for each fold the same number of features chosen by PFE. For each fold, each algorithm is used to select a subset of the original features.

In order to generate reproducible results, all algorithms that exploit pseudorandom elements in their training process have been set with a fixed seed. Unless differently specified, each algorithm uses its default parameters as defined in [START_REF] Pedregosa | Scikit-learn: Machine learning in python[END_REF][START_REF] Li | Feature selection: A data perspective[END_REF]. Each dataset has been standardized removing the mean and scaling to unit variance (StandardScaler [START_REF] Zill | Advanced engineering mathematics[END_REF]).

Redundancy detection

The MADELON dataset proposed in [START_REF] Guyon | Design of experiments of the nips 2003 variable selection benchmark[END_REF] was specifically designed to challenge feature selection algorithms in detecting redundant features. Features generated by MADELON can be informative (d i ), repeated (d r ), or redundant (d c ). The algorithm generates clusters of points normally distributed about vertices of an hypercube in a subspace of dimension d i and assigns an equal number of clusters to each class. Then it stacks d c linear combinations of the informative features followed by d r duplicates, drawn randomly with replacement from the informative and redundant features. All the remaining features (d n ) are random noise

(d n = d -d i -d c -d r )
. This benchmark dataset is used to assess the ability of PFE in detecting redundant features. For this experiment the number of informative features is set to 150 as well as the number of redundant and repeated features (d i = d c = d s = 150). The total number of features is set to d = 500, thus d n = 50 features are just random noise. The task is to detect informative, redundant, duplicate, and noisy features in order to correctly classify clusters of samples on hypercube vertices. Experimental results are shown in Figure 1. Once feature selection algorithms are fitted on a training fold, an instance of RidgeClassifier and of DecisionTreeClassifier are used to assess the quality of the selection on the test set. The resulting F 1 score [START_REF] Van Rijsbergen | Information retrieval. 2nd. newton[END_REF] is then compared to the one obtained by training on the same fold but using all the original features. In this way, the performance of all techniques can be evaluated with respect to a fair baseline (see Figure 1, top). Except for lap score, PFE resulted as the fastest approach. The most interesting result of the simulation is represented by boxplots in Figure 1, bottom. They show for each kind of feature (informative, redundant, duplicate, and noisy) the percentage of features retained by each algorithm. Notably, PFE and MCFS preserve most of the informative features while discarding most of the redundancy. However, MCFS is the worst technique in terms of duplicate detection, whereas PFE is the best one together with NDFS and lap score. It should be remarked, anyway, that PFE retains all noisy features. Yet, it is not surprising at all as the approach is not designed to get rid of random noise. In authors view, PFE is not meant to be used alone but combined with other complementary feature selection approaches.

Cross-task benchmarks

The ability of feature selection algorithms in tackling different kind of machine learning problems is assessed using four benchmark datasets taken from the OpenML [START_REF] Vanschoren | Openml: Networked science in machine learning[END_REF]. Table 1 highlights the main characteristics of the four datasets. The first two (gas-drift and isolet) are used to test classification and clustering performances while the latters are employed for regression. The quality of the selections is assessed over a 10-fold cross-validation. Once fitted on a training fold, each algorithm provides a selection of the original features to an instance of a machine learning model on the filtered training set only. RidgeClassifier and DecisionTreeClassifier are used to evaluate classification performance in terms of F 1 score. AgglomerativeClustering [START_REF] Ward | Hierarchical grouping to optimize an objective function[END_REF] and KMeans [START_REF] Steinhaus | Sur la division des corp materiels en parties[END_REF] are employed to assess clustering performances through the Silhouette coefficient [START_REF] Rousseeuw | Silhouettes: a graphical aid to the interpretation and validation of cluster analysis[END_REF]. The number of centroids for k-means is chosen as twice as much as the number of classes. For regression Ridge and DecisionTreeRegressor are used to measure the coefficient of determination (R2 score, [START_REF] Steel | Principles and procedures of statistics[END_REF]). Once collected, performance scores are compared to the ones obtained by training on the same folds but using all the original features. In this way, the performance of all techniques can be evaluated with respect to a fair baseline. Figures 2,3, and 4 show the results in terms of performance metrics and training time. As mentioned before, all the other feature selection algorithms are set in order to provide for each fold the same number of features chosen by PFE, thus yielding a fair comparison.

Compared to state-of-the-art techniques, PFE is among the fastest solutions together with RFE and lap score. Notably, RFE is not used for clustering as it is a supervised algorithm, thus it cannot be employed for unsupervised tasks. Despite its unsupervised nature, PFE often matches RFE performances and sometimes provides even better solutions (i.e., Mercedes dataset using DecisionTreeRegressor). The efficiency of PFE with respect to other unsupervised approaches is revealed on the largest dataset (gas-drift) where it is faster by a few order of magnitudes. Moreover, even when PFE performances appear to be slightly worse than others (i.e., Mercedes using Ridge), it may be sufficient to change the downstream predictor (i.e., the performance looks much better when DecisionTreeRegressor is used). Indeed, by construction, PFE performs feature selection such that the information loss is almost negligible.

Conclusions

In this paper, a novel feature selection approach named Predictable Features Elimination has been introduced. At the heart of the methodology lies the idea that features whose value can be easily predicted based on the values of other features of the same sample, probably contain mostly redundant information. The algorithm iteratively trains a machine learning model using one of the features as a target, and if the quality of the model is above a user-defined threshold (e.g., V > 0.9), the feature is removed. After all features have been treated in this way, the remaining ones will constitute the final feature set. Not relying on target labels or values, PFE can be used for supervised, unsupervised, or semi-supervised machine learning problems.

Experimental results prove PFE to be competitive with state-of-the-art feature selection algorithms, on a set of non-trivial classification, regression, and clustering benchmarks. The main drawback of the approach is the impossibility of removing uninformative, but hard-to-predict features, for example those including completely random values: In most cases, however, such features are filtered out by subsequent machine learning algorithms applied to the data, as they have no correlation with the objectives. Future works will investigate the performance on PFE on a wider range of benchmarks, and explore the possibility of using a similar idea on samples, to uncover coresets and potentially perform dataset compression. Finally, particular focus

Fig. 1 .

 1 Fig. 1. Classification performances on the MADELON dataset (top two); feature selected for the MADELON dataset (bottom four).

Fig. 2 .

 2 Fig. 2. Classification results on benchmark datasets.

Fig. 3 .

 3 Fig. 3. Regression results on benchmark datasets.

Fig. 4 .

 4 Fig. 4. Clustering results on benchmark datasets.

Table 1 .

 1 Benchmark datasets.

	dataset samples features classes reference
	gas-drift 13,910	128	6	[27]
	isolet	7,797	617	26	[8]
	Mercedes 4,209	377	-	[7]
	crime	1,994	127	-	[24]

https://github.com/glubbdubdrib/predictable-feature-elimination


Eliana Giovannitti

Sayyidshahab Nabavi

Giovanni Squillero

Alberto Tonda

Exploiting Artificial Swarms for the Virtual Measurement of Backlash in Industrial Robots

Keywords: backlash, robotic manipulator, real-world, industry application, swarms-based approaches

The backlash is a lost motion in a mechanism created by gaps between its parts. It causes vibrations that increase over time and negatively affect accuracy and performance. The quickest and most precise way to measure the backlash is to use specific sensors, that have to be added to the standard equipment of the robot. However, this solution is little used in practice because raises the manufacturing costs. An alternative solution can be to exploit a virtual sensor, i.e., the information about phenomena that are not directly measured is reconstructed by signals from sensors used for other measurements.

This work evaluates the use of bio-inspired swarm algorithms as the processing core of a virtual sensor for the backlash of a robotic joint. Swarm-based approaches, with their relatively modest occupation of memory and low computational load, could be ideal candidates to solve the problem. In this paper, we exploit four state-of-the-art swarm-based optimization algorithms: the Dragonfly Algorithm, the Ant Lion Optimizer, the Grasshopper Optimization Algorithm, and the Grey Wolf Optimizer. The four candidate algorithms are compared on 20 different datasets covering a range of backlash values that reflect an industrial case scenario. Numerical results indicate that, unfortunately, none of the algorithms considered provides satisfactory solutions for the problem analyzed. Therefore, even if promising, these algorithms cannot represent the final choice for the problem of interest.

I. INTRODUCTION

Industrial robots are used by manufacturing companies to perform repetitive, difficult or even dangerous tasks [START_REF] Edwards | Robots in industry: An overview[END_REF]. Robotics manipulators are made of a sequence of rigid links, connected by movable joints. These joints are actuated by electrical motors and include a set of gears that are designed to deliver power from the motor to the link as efficiently as possible. In engineering, Backlash is the name given to the lost motion in a mechanism caused by gaps between the parts [START_REF] Stein | Estimation of gear backlash: Theory and simulation[END_REF]; in industrial manipulators, that is, in the context of gears and gear trains, the backlash is the amount of clearance between mated gear teeth (Fig. 2). Backlash is dependent on Authors are listed in alphabetical order. wear and increases over time [START_REF] Yang | Gear backlash detection and evaluation based on current characteristic extraction and selection[END_REF]. It causes vibrations to arise and eventually impair the robot positioning accuracy degrading the performance. The immediate effect is a reduction in the quality of the pieces produced by the robot, and in the long term the breakage of the joint can occur. To prevent this scenario industries target maintenance interventions with great care. On the one hand, they try to avoid unnecessary components substitutions and limit the number of stops; on the other hand, they don't want to come up with the breakdown of a joint that may cause as a stop in the whole production line bringing significant economic losses. Therefore, to guarantee the quality of the production and to predict failures, an assessment of the value of the backlash is fundamental. A direct measure of the backlash can be obtained using two coupled encoders, one on the motor side and one on the load side, or using an accelerometer. Unfortunately, these sensors are not part of the standard equipment of industrial robots, and adding them to the robot means higher costs of the final product. Moreover, it can be difficult to find free space inside the robotic joint to place the extra sensors. A solution for these problems can be the use of virtual sensors. A virtual sensor is a piece of code which serves as physical sensor. It uses the information available from other onboard measurements, and infer an estimate of the quantity of interest 978-1-7281-8393-0/21/$31.00 ©2021 IEEE that can not be directly measured. The measure is achieved relying on mathematical models and relationships between different variables. Fig. 2: Backlash in mating gears.

In [START_REF] Giovannitti | Virtual measurement of the backlash gap in industrial manipulators[END_REF], authors try to estimate the value of the actual measurement of the joint backlash by using measures from the onboard motor encoder. The measure of the motor position is carefully analyzed looking for information about the backlash status. In particular, the analog signal is scanned to find the match with a peculiar disturbance pattern related to the backlash presence, see Fig. 3. This pattern is obtained and characterized through simulation, and is taken as the signature of the backlash in the encoder signal. A generic formula ( 1) is developed to model it, and an estimate of the underlying backlash is obtained by fitting the formula against the measured data. The measure of the backlash is simply given by the measure of the amplitude of the oscillation, but the fit of the whole formula is required to discern the specific disturbance caused by backlash from any other disturbance source. In more Fig. 3: Backlash disturbance on the reference signal (blue) and how its appearance changes at increasing backlash gap. details, the formula of the expected disturbance on the encoder signal is

d b (t) = A e -(t-t1)τ sin ω(t) t 1 < t < t 2 0 otherwise ( 1
)
where A and τ are the parameter related to the amplitude decay; t 1 is the starting time of the oscillation, and t 2 is the ending time, that is,

d b (t) = 0 if t ∈ [t 1 , t 2 ].
In [START_REF] Giovannitti | Virtual measurement of the backlash gap in industrial manipulators[END_REF] the optimization phase is carried out by an Evolutionary Algorithm that provides accurate and repeatable results. Unfortunately, the algorithm has proved to be time consuming, and demanding in terms of memory occupation and computational load. From this emerged the idea of this paper: try to use easier and lighter to implement optimization algorithms.

In this work, we analyze the use of different swarm-based meta-heuristics for the core of the virtual sensor, that is the fitting of equation ( 1) with the real data, and evaluate their performance in terms of precision and accuracy.

The paper is organized as follows: Section II describes the four algorithms analyzed and the bio-inspired behavior models at the basis of their functioning; Section III provides an overview of the backlash estimation problem and how it has been formalized in terms of an optimization problem; Section IV is about the experiments performed and the results discussion. Finally, the last section contains the conclusions on the work carried out.

II. SWARM-BASED APPROACHES

This paper is focused on four state-of-the-art nature-inspired algorithms, all emulating typical behaviours of the animal world. They are: the Dragonfly Algorithm, the AntLion Optimizer, the Grasshopper Optimization Algorithm, and the Grey Wolf Optimizer. The algorithms were all developed by the same author, Seyedali Mirjalili, and are freely downloadable from the Web. The logic behind them is described as follows.

A. Ant Lion Optimizer

Ant Lion Optimizer (ALO) was developed in 2015, [START_REF] Mirjalili | The ant lion optimizer[END_REF]. The idea was inspired by the hunting mechanism of the antlion (a species of insect in the family Myrmeleontidae). The hunting strategy of this animal is based on the use of traps it creates to capture its victims: the ants. The antlion digs a cone-shaped pit in the sand, then conceals under the trap and waits for the victim. As soon as the insect finds the prey in the trap, it starts throwing sand towards the edge of the trap making the victim slump down to the bottom of the trap. After having its meal, it prepares the trap for the next time.

The ALO algorithm takes place in the following 4 main steps: the definition of a random position and a random movement for the ants, the establishing of traps, the enmeshment of ants in the trap, the capture of the victims. Ants and antlions are associated to the possible solutions, at every iteration their fitness is evaluated. Ant move randomly, while antlions are more stationary, hidden in their traps, in defined points of the search space. The position of an antlion represents the best solution found so far in the area around it. This area correspond to the size of the antlion trap, and is proportional to the fitness of the antlion: the greater the fitness, the greater the area of the trap. Each ant is associated to a corresponding antlion, and can only move inside the area defined by the trap. The ant's task is to explore this area looking for a solution with a fitting higher than the fitting of the antlion. If this happens, the position of the antlion is updated with the position of the ant and the radius of the trap is increased accordingly. This phase correspond to the capture of the prey. Finally, the antlion with the best fit is considered as an elite, and when optimization ends, the position of this antlion represents the solution of the problem.

In ALO, the random walk followed by ants is modeled by

x(t) = [0, S(2r(t 1 ) -1), ..., S(2r(T ) -1)] (2)
where S(•) is the cumulative sum, t i counts the iterations, T indicates the maximum number of iterations of the algorithm, and r(t) is a stochastic function that returns 0 or 1 value with equal probability. The position of the ants during optimization is stored in the following matrix:

M A =      A 11 A 12 . . . A 1d A 21 A 22 . . . A 2d . . . . . . . . . . . . A n1 A n2 . . . A nd      (3)
where A ij is the position in the j-th dimension of the i-th ant. The number of the ants is indicated by n, while d is the problem dimension. Given f (•) as the fitness function, the fitness values for each ant during optimization are stored in the matrix

F A =      f ([A 11 , A 12 , ..., A 1d ]) f ([A 21 , A 22 , ..., A 2d ]) . . . f ([A n1 , A n2 , ..., A nd ])      (4)
These two matrices are updated at each step of the optimization. Similar matrices, M AL and F AL , are defined for antlions. To set the position for the antlions a roulette wheel selection is applied. Antlions are selected based on their fitness values, a higher chance is given to antlions with higher fitness value. To model the falling of the prey in the cone shaped trap, the radius of random walk used in each step is decreased following these equations

C t = C t I d t = d t I (5)
where, C t and d t indicates the minimum and the maximum, respectively, of all variables at iteration t, I is the decreasing ratio calculated as I = 10 w t T . In this last equation, T is the maximum number of iteration and w is a constant defined based on the current iteration. The value of w gradually increases as iterations proceed, controlling the accuracy level of the exploitation stage.

To model the effect of antlion's ambush on the random movement of the ants the following equations are used

C t i = Antlion t j + C t d t i = Antlion t j + d t (6)
in which Antlion t j is the position of antlion j at iteration t. Moreover, the random movement of an ant around the corresponding antlion is also affected by the position of the elite; so the position of the i-th ant at the t-th iteration, Ant t i , is expressed as

Ant t i = R t A + R t E 2 (7)
where R t A and R t E are the random walks around the antlion and the elite considered at the same iteration.

B. Grasshopper Optimization Algorithm

The Grasshopper Optimization Algorithm (GOA) [START_REF] Mirjalili | Grasshopper optimization algorithm for multi-objective optimization problems[END_REF] is a swarm-based algorithm that imitates social behavior of grasshoppers seeking for food.

In their lifetime, these insects pass through two main phases: nymph and adulthood. While nymph, they have no wings so they move slowly and eat all the vegetation on the floor. When adult, they develop wings so they can move fast and fly covering long distances to find a new place with food. The swarming behavior characterizes both nymph and adulthood. The two modality grasshoppers use to search for food can be easily compared to the exploitation and exploration phases in an optimization algorithm.

In GOA, grasshoppers represent the possible solutions. The insects move in the search space, and the fittest one at each step of the optimization is elected as the food position, i.e. the best solution so far. The food position is passed to the other individuals of the swarm requiring them to move toward it. While other grasshoppers move, they can find a fitter position than the previous one. This last position then becomes the position for the food in the new iteration. While the optimization progresses, the algorithm continuously finds better solutions and converges to the optimum.

In GOA the position, X, of i-th grasshopper is given by the formula

X i = S i + G i + A i (8)
where S i , G i , A i represent the three forces that act on the insect behavior. The social interaction, the gravitational attraction, and the wind flow, respectively. The social interaction term, which has the highest impact on the movement, is expressed as

S i = N j=1,j =i s(d ij ) d ij . (9)
In the formula, d ij is the distance between grasshopper i and grasshopper j. It is computed as

d ij = |x j -x i |. While d ij = xj -xi
dij is the unit vector from the i-th grasshopper to j-th one. The function s(.) represents the power of social interaction, i.e, attraction or repulsion, among grasshoppers and is expressed as

s(r) = f e -r l -e -r (10)
where f represents attraction intensity and l is the the scale of attractive length. Thanks to this formula, the space between two grasshoppers is divided into a repulsion zone (the inner of the hypersphere centered in the i-th grasshopper), an attraction zone (the outer of the hypersphere), and a comfort zone (the surface of the hypersphere) where there is neither attraction nor repulsion. Changing the value of l and f changes the social behavior of grasshoppers. The gravitational term in equation ( 8) is given by

G i = -g e g ( 11
)
where g is the gravitational constant, and e g is the unity vector in the direction of the earth center.

The last component in the position formula ( 8) is the wind drift, A i , computed as

A i = u e w ( 12
)
where u is a constant drift, and e w is a normalized vector in wind direction. From all the above, the equation that defines the position of the i-th grasshopper can be rewritten as

X i = N j=1,j =i s(|X j -X i |) x j -x i d ij -g e g + u e w ( 13
)
where N is the total number of grasshoppers. However, this formula is a good model for the behavior of the elements of the swarm but needs a further modification to be used for an optimization algorithm. Indeed, the formula quickly leads, and then holds, the grasshoppers in the equilibrium position, without making them converge in a single point. So, it is modified as follow

X d i = c N j=1,j =i c ub d -lb d 2 s(|x d j -x d i |) x j -x i d ij + T d (14)
where c is a coefficient which decreases as the iterations proceed shrinking the repulsion region between grasshoppers, ub d and lb d are the upper and lower bound in the d-th dimension, T d is the best value found so far in dimension d, and s(.) is approximately the same function used in [START_REF] Mirjalili | Dragonfly algorithm: A new meta-heuristic optimization technique for solving single-objective,discrete, and multi-objective problems[END_REF].

The G i contribution is not considered, and it is presumed that the wind (i.e., the A i term) is always in the direction of T d . The c factor appears twice in the formula because is also used to act on the exploration and the exploitation behaviour. While decreasing, c reduces the movements of grasshoppers around the target, promoting the exploitation at the expense of exploitation. It is updated with the following equation

c = c max -l c max -c min L (15)
where c max is the maximum value, c min is the minimum value, l indicates the current iteration, and L is the maximum number of iterations.

C. Grey Wolf Optimizer

The Grey Wolf Optimizer (GWO) is inspired by the social and hunting behavior of grey wolves in nature, [START_REF] Mirjalili | Grey wolf optimizer[END_REF].

Grey wolves live in packs and respect a strict social hierarchy which divides them into 4 main categories: alpha, beta, delta, and omega. The dominant wolf is the alpha, it leads the pack. The beta and delta wolves help the alpha in making decisions about pack activities or accomplice main tasks like hunting or watching the boundaries of the territory. All other wolves in the group are omega individuals. The hunting strategy of grey wolves can be summarized in three phases: searching for prey, encircling prey, and attack prey.

The GWO logic mimic these behaviours. Each wolf is associated to a possible solution. At every iteration the solutions are evaluated and fittest one is taken as alpha wolf, while the second and third best solutions became beta and delta wolves. These three solutions are used to define an estimate of the position of the prey, and this information is then used to update the position of wolves for the next iteration. Wolves have to converge to the prey while diverging form each other. The exploration (diverging) or the exploitation (converging) behavior of wolves can be mitigated or enhanced acting on specific coefficients in the position formula.

The generic position of a wolf around a prey can be expressed by using two equations:

→ D = | → C. → X p (t) - → X(t)| (16) → X(t + 1) = | → X p (t) - → A. → D|. (17)
In these equations, t is the current iteration, → X p indicates the position vector of the prey, while → X is the position vector of a grey wolf.

→

A and → C are coefficient vectors given by the formulas

→ A = 2.a. → r 1 ( 18
) → C = 2. → r 2 (19)
where → r 1 , and → r 2 are random vectors in [0,1], and a is a coefficient that is linearly reduced from 2 to 0 while the optimization proceeds. Since the position of omega wolves is updated with respect to the alpha, beta, and delta ones, the following equations are used

→ D α = | → C 1 . → X α - → X| (20)
→ D β = | → C 2 . → X β - → X| ( 21
) → D δ = | → C 3 . → X δ - → X|. ( 22
) → X α , → X β , →
X δ indicates the position of alpha, beta and delta wolves respectively, and C 1 , C 2 , C 3 are random vectors. To estimate the approximate distance between the current solution and alfa, beta and delta wolves, equations (20), (17), and (22) are applied as follow:

→ X 1 = → X α - → A 1 . → D α (23) → X 2 = → X β - → A 2 . → D β ( 24
) → X 3 = → X δ - → A 3 . → D δ (25)
Then, the position of the current solution can be updated for the next iteration, and becomes

→ X(t + 1) = → X 1 + → X 2 + → X 3 3 . ( 26
)
It can be observed that two vectors,

→

A and → C, are defined. These two vectors are random and are used to tune the strength of the exploration and exploitation phases in the GWO algorithm.

D. Dragonfly Algorithm

The last algorithm evaluated in the present work is the Dragonfly Algorithm (DA), [START_REF] Mirjalili | Dragonfly algorithm: A new meta-heuristic optimization technique for solving single-objective,discrete, and multi-objective problems[END_REF]. The DA is an optimization method inspired by the swarming behavior of dragonflies. They exhibit two different attitudes while flying: a static feeding swarming, and a dynamic migratory swarming. These two habits loosely mimic the characteristics of the two phases of an optimization process: exploration and exploitation.

The base principles used in this algorithm are the same swarm principles identified by Reynolds in [START_REF] Raynolds | Flocks, herds, and schools: A distributed behavioral model[END_REF] and [START_REF] Raynolds | Steering behaviour for autonomous characters[END_REF]. Elements in a swarm follow three fundamental rules while moving:

• Separation, which means avoid collisions between individuals; • Alignment, which indicate individuals should have the same velocity; • Cohesion, which refers to the tendency of individuals to converge to the center of the mass of the swarm. Two more principles can be added to the previous three:

• Attraction towards food, which is the main aim of the swarm; • Distraction outwards enemy. In the following, a model for each principle is provided. The separation term is given by

S i = - N j=1 X i -X j (27)
where the position of the current individual is indicated by X i , N is the number of agents in X i 's neighborhood, and X j is the position of the j-th neighboring individual. A radius r is used to define a region of visibility for the considered dragonfly, all dragonflies inside this region are taken in account in the calculation of S i . The same visibility rule is used for the other behavioural terms. The alignment is computed as

A i = N j=1 V j N ( 28
)
where V j represents the velocity of the j-th neighboring individual. While the cohesion is given by

C i = N j=1 X j N -X i . ( 29
)
Attraction towards the food source and distraction outwards the enemy are computed as

F i = X + -X i ( 30
)
E i = X -+ X i (31)
where X + indicates the position of the food source, while X - is the position of the enemy. Agents move within the search space and their position is updated at every iteration with the rules:

X t+1 = X t + ∆X t+1 (32)
∆X t+1 = (sS i + aA i + cC i + f F i + eE i ) + w∆X t (33)
The vector ∆X is the velocity vector of the dragonfly and contains information about the direction and speed of the movement. The constant w is an inertia weight, and t is the iteration counter. The coefficients s, a, c, f , e are used to weight the different behaviours. In the first steps of the optimization process the dragonflies fly in a dynamic swarm. The fly is highly coordinated, dragonflies speeds are aligned and cohesion and separation are high. This is the best organization for exploring the search space. While the optimization proceeds, the swarm behavior changes and gradually becomes a static swarming. In this last stage the cohesion increases while alignment and separation decrease, so that the dragonflies can converge towards the prey exploiting the search space. The transition between the two flight modes occurs by acting on the radius r and on the coefficients s, a, c, f, e. Their values are modified with the increasing number of iterations. A last note should be added about food and enemy position in DA.

They are updated at the beginning of every iteration, when the fitting of all dragonfly is evaluated. The position of the fittest one is taken as the food location, while the position of the worst one is considered as the place where the enemy is.

III. METHODOLOGY

The simple logic and the easy equations of the presented algorithms have made them of interest for the industrial problem addressed. Since the idea is to implement a virtual sensor, that is a software sensor, the occupation of memory and the computational effort are taken carefully into consideration. These aspects are even more important when dealing with industrial devices where software resources are often tailored to the primary needs of the system. In order to use swarm algorithms to address the backlash problem this last must be translated into an optimization problem. Details on this are presented below.

The backlash phenomena arises inside the joint, but its effects propagate both to the link and to the motor that are connected at the input and output ends of the joint, see Fig. 4. In standard industrial robots there are no sensors in the joint, nor on the link. While an encoder is always present on the motor. The encoder is used in position and speed control loops of the motor itself. By detecting the consequences on the motor of the backlash in the joint, is possible to derive information on the backlash presence in the joint itself. However, the detection is not so simple because these effects can reach the motor strongly attenuated and be added or confused with other sources of disturbance.

The main idea of this virtual sensor is to identify a model as the signature for the backlash effect on the motor and look for this pattern inside the signal from the encoder. When the model is detected, its parameters are measured end an estimate of the backlash level in the joint is derived. The greater the amplitude of the disturbance, the greater is the backlash.

The matching between the model and the signal, and the measurement of the model parameters are performed through the optimization algorithm. All the information about the cost function, the parameters, the search space and the constraints of the problem are described in what follows.

A. The objective function

For the problem of interest, Root-mean-square error (RMSE) between data and model is the objective function to minimize. Since data show a repetitive pattern, made of positive and negative sequential oscillations, the disturbance signature (1) has been slightly modified to replicate this behavior too. resulting formula is

h(t, A, t 0 , τ, ω, v t , t 1 , t 2 , T, T w ) = v t + 12 i=1 f (t, A, t 0 + i • T w , τ, ω, t 1 , t 2 , T ) (34)
where v t is the commanded motor speed, and the function f is the sequence of a positive and a negative, T shifted, disturbance oscillation

f (t) = d b (t -t 0 ) -d b (t -t 0 + T ). (35)
The t 0 parameter is the disturbance starting time while T w is the time period for the f (t) function. So, it turns out that the objective function is

RM SE = Σ N i=1 v(t) -h(t, A, t 0 , τ, ω, v t , t 1 , t 2 , T, T w ) 2 N
where N is the number of samples in the dataset.

B. The parameters

The parameters to be identified are

X = [A, t 0 , τ, ω, v t , t 1 , t 2 , T, T w ] (36)
The vector dimension is 9 but, since T and T w have fixed and known in advance values, the problem dimension becomes 7. Considering that all the parameters have a physical meaning and that some information about the system under test is apriori available, it is possible to use it to define upper and lower bounds for the parameters. In particular, boundaries for amplitude related parameters of the model, A and c, are proportional to the amplitude of the signal to be analyzed. We know in advance they have increasing amplitude because the disturbance superimposed on it is growing. Furthermore, these parameters must have a definite relationship to one another to ensure that the aspect of the function respects the expected one. In particular the constraints are defined for t 0 , t 1 , and t 2 and are: 0 < t 0 < (max(t) -t 1 ), t 2 < t 1 < T and 0 < t 2 < pi/(2ω).

At the end of every iteration, when the position of the agents is updated for the next iteration and after the check to verify if the new values fall within the defined domains, the constraint check is also performed. If the value of a parameter is outside the interval defined by the constraints, then the value is saturated to the extreme of the allowable range.

To perform a comparison, the algorithms are run under the same conditions: same number of agents and iterations, same constraints and same input data. All algorithms used are applied in their original Matlab implementation. No changes are made to the code other than minimal customization to fit the problem.

C. The datasets

Each algorithm is run on the same test set. The test set comprises 20 signals, all affected by a different level of backlash.

Signals are obtained with the help of simulation. The Matlab/Simulink environment is used to model the robotic joint and record the motor encoder signal when the system is affected by an increasing value backlash (see Fig. 3). Dataset1 contains the signal that corresponds to the lowest level of backlash, while Dataset20 contains the data corresponding to the highest backlash. Fig. 5 shows the datasets and the corresponding values for the backlash (in blue) and for the disturbance amplitude (in red). The datasets correspond to the 20 backlash values that span the interval [δ min ; δ max ] = [0.0001; 0.004]radians with the constant step ∆δ = 0.0002radians. Since it represents the minimum distance between two possible measures, the step was also taken as the accuracy value required for the virtual sensor.

D. The optimization

The number of agents and the number of iterations used for the four algorithms is the same. The best setting is identified in 300 agents and 200 iterations. The value is chosen taking the DA algorithm as a reference and following a trial and error strategy. Since the aim of our study is to compare the performances of the different optimization methods, it is decided to proceed under the same conditions for all the algorithms and use [300,200] as the common setting.

The workstation used for carrying out the result is equipped with a 2.70 GHz processor, 16 GB of RAM, and the system used is Microsoft7. To run the proposed algorithm the Matlab software is used.

IV. EXPERIMENTS

To evaluate the performance of the different algorithms with respect to the problem of interest, an extensive test campaign is conducted. The four algorithms are run on the backlash measure problem described in Section I, and the identification is performed on the datasets described in Section III. Since preliminary tests have shown that in the first two datasets the backlash disturbance is too small and cannot be properly identified, these dataset are ignored in the analysis. The remaining 18 datasets have served as reference for the assessment. Thirty independent experiments are performed for each set while all conditions are kept as constants. The same setting is used for all tests and for all the algorithms. Furthermore, to provide each algorithm with equal opportunities, the starting solution is always randomly chosen.

A collection of all results is shown in Fig. 6. Graphs in the upper part summarize the values of the backlash disturbance amplitude provided by the different optimization methods. Dots stacked in each column correspond to the 30 outcomes of experiments on a same dataset. The red line in the graph connects the expected values for the disturbance amplitude. While the red dotted lines represent the maximum allowable error for the measure, i.e. the accuracy on δ reported on the value of A. Since the relationship between A and δ is not linear, the accuracy on A is taken as the minimum distance between two consecutive values of A plotted in Fig. 5. Looking at the figures, it can be noted that the distribution of results is always widening while moving towards larger values of the backlash, i.e., going left to right in the plots. The reason behind this is not the increased difficulty of the problem, but is the way the problem search space is defined. In fact, amplitude related parameters of the model, A and c, have boundaries that grow while moving through the datasets on the right leading to an increase in their search space, see Section III-B. It is also possible to note that it is more likely that the swarm optimization algorithms tend to underestimate, rather than overestimate, the solution. It can be attributed to the fact that the fitting result is sometimes wrong. When fitting fails, only the second and smaller half of the decreasing oscillation is fitted, while the first and largest part is lost. In such cases an underestimated solution is provided.

A different perspective on the same data is given by the box and whisker diagrams, in Fig. 7. They present in quantitative terms the dispersion of the results around the expected value. Diagrams give clear information about GOA, ALO and DA performance, but unfortunately are not significant for GWO. This is because the particular shape of the dispersion of GWO results, with the majority of them lying very close to the red expected value and the wide empty band separating the remaining data (see GWO graph in Fig. 6), makes the box plot misleading. So, even if box charts show ALO as the best in terms of precision, a more in-depth analysis returns GWO as the best of all.

Since comparing the outcome of stochastic algorithms is not trivial, as just examining the mean performance is not enough, the results are also evaluated through a statistical test procedure. The series of 30 results obtained from experiments on the single dataset is treated as a set of samples from a continuous distribution. The two-sample Kolmogorov-Smirnov test [START_REF] Smirnov | On the estimation of the discrepancy between empirical curves of distribution for two independent samples[END_REF] is used to check the hypothesis that two different series of results come from the same distribution. In particular, the test is run with α = 0.05 as the threshold for significance of distinctiveness. The results are shown in Table II. The reported value is the mean of the fitting error on the 30 repeated experiments. Italic is used for results that are likely form the same distribution, and therefore not significant for the comparison. For the majority of cases, however, it is possible to accept the hypothesis that the data collected come from different distributions. In these cases the best result is indicated in bold. Again, GWO appears as the best of the four algorithms. For all the above, it can be stated that the GWO algorithm has the best performances with respect to the backlash problem. While, DA is the worst. Nevertheless, none of the algorithms analyzed presents performances that can meet the requirements of precision and repeatability of our sensor. Fig. 6 clearly shows the solutions lying outside the allowable band delimited by the desired accuracy. Therefore, at the present stage of development, these algorithms cannot be used for the implementation of the virtual sensor for the backlash.

V. CONCLUSIONS

The paper presented an application of four bio-inspired swarm algorithms to a real industrial problem. A extensive test campaign was conducted and a thorough comparison of results and performances of the algorithms was performed. Key aspects as accuracy, precision, and ease of implementation were taken into account, focusing on the possible use of the swarm algorithms for the implementation of a virtual sensor. The problem of interest, i.e., the matching between a disturbance signature and a signal, has proved difficult to solve, and the considered algorithms were found to be unsuitable to find its solution. The presence of many local minima has been observed, and none of the algorithms showed sufficient exploration capabilities not to get stuck in one of the many local solutions. It has been verified that, even showing some advantages in terms of memory occupation and ease of implementation, the bio-inspired swarm algorithms can not be used for the backlash sensor implementation. They were not able to guarantee the level of reliability and accuracy required by the sensor.

Fig. 1 :

 1 Fig. 1: A robotic manipulator (Comau SMART NJ 650).

Fig. 4 :

 4 Fig. 4: Simplified model of a robotic joint.

Fig. 5 :

 5 Fig. 5: Backlash and expected disturbance amplitude on the considered datasets.

Fig. 6 :

 6 Fig. 6: Collection of all the results of the experiments.

Fig. 7 :

 7 Fig. 7: Box and whisker diagram.

TABLE I :

 I Parameters variability range

	Symbol		Min value		Max value		Units
	A	-1.3 * (	max v(t)-min v(t) 2	) 1.3 * (	max v(t)-min v(t) 2	)	rpm
	t 0		min t		max t		s
	τ		5		30		-
	ω		2π		2π • 40		rad/s
	vt		min v(t)		max v(t)		rpm
	t 1		0		0.204		s
	t 2		0		0.204 2		s
	T		0.204		0.204		s
	Tw		0.408		0.408		s

TABLE II :

 II Kolmogorov-Smirnov test results

	Datasets	ALO	GOA	DA	GWO
	3	2.777	4.781	3.227	2.637
	4	3.331	6.169	4.558	3.278
	5	4.560	7.214	6.039	3.327
	6	5.752	9.511	7.422	3.630
	7	6.728	8.946	9.277	4.250
	8	6.622 11.628	9.363	5.400
	9	7.580 12.163 11.553	5.817
	10	8.816 13.065 11.412	7.310
	11	9.456 14.159 13.399	7.957
	12	8.647 16.212 12.846	8.268
	13	10.218 16.571 14.409	8.844
	14	10.729 18.311 14.321	9.289
	15	10.957 17.419 16.261 10.411
	16	10.427	19.766 16.709 10.509
	17	11.301 20.982 16.921 10.521
	18	11.932 20.055 18.920 11.017
	19	13.628	20.69 17.799 12.576
	20	13.660 22.126 19.150 11.214