Alice Nicolaï
email: alice.nicolai@cmla.ens-cachan.fr

Julien Audiffren

Model-Space Regularization and Fully Interpretable Algorithms for Postural Control Quantification

As falls prevalence increases with the aging of the population, early detection of balance degradation is of great importance for efficient prevention and treatment. This work addresses the problem of quantifiying static balance with fully interpretable learning algorithms. Our approach relies on a heuristic based variant of the aggregation of weak classifiers constrained with a new model-space regularization combined with a family of interpretable features. In our experiments, these models outperforms their regular alternative, opening promising new research directions.

Introduction

In the recent years, Machine Learning methods and algorithms have been successfully used in a wide range of problems, ranging from robotics [START_REF] Kober | Reinforcement learning in robotics: A survey[END_REF] to recommander systems [START_REF] Goodfellow | Deep learning[END_REF]. In particular, possible applications in Healthcare and Medical Research have gained significant interest recently, as the ageing of the populations will increase the need for efficient and widely available healthcare and treatments [START_REF] Christensen | Ageing populations: the challenges ahead[END_REF]. Notably, the screening and early diagnosis of diseases and conditions is one of the focal point of attention for Machine Learning methods, and Deep Learning approaches have already lead to breakthtough results in Medical Imagery (see e.g. [START_REF] Sirinukunwattana | Locality sensitive deep learning for detection and classification of nuclei in routine colon cancer histology images[END_REF]). Other methods have been proven valuable for computer-aided diagnosis [START_REF] Li | Improve computer-aided diagnosis with machine learning techniques using undiagnosed samples[END_REF], in particular for conditions where the early symptoms may be difficult to assess.

One domain in which early diagnosis is key is balance deterioration. Indeed, such impairments strongly increase the risk of falling, which can have devastating consequences -including serious injuries, dependency and death -particularly for the elderly population [START_REF] Sattin | The incidence of fall injury events among the elderly in a defined population[END_REF]. Consequently, fall prevention in the elderly population is of paramount importance, particularly as each year more than a third of population 65 years-old and older faces a fall [START_REF] Tinetti | Preventing falls in elderly persons[END_REF]. A key ingredient to fall prevention is the early detection of balance deterioration, which can then be treated by prevention strategies such as reeducation [START_REF] Van Diest | Exergaming for balance training of elderly: state of the art and future developments[END_REF]. Accurate quantification of balance control can be of significant assistance to such detection, and can also be used to evaluate the performances of prevention and recovery strategies.

Quantifying Balance. Balance is achieved by a complex synergy of many biological systems, including the visual, proprioceptive, vestibular and central nervous system [START_REF] Peterka | Role of somatosensory and vestibular cues in attenuating visually induced human postural sway[END_REF] -and therefore is difficult to thoroughly evaluate. Generally, balance is considered as the combination of a static and a dynamic component [START_REF] Perrin | Influence of visual control, conduction, and central integration on static and dynamic balance in healthy older adults[END_REF]. Static balance is defined as the capacity to stand still without falling whereas dynamic equilibrium is the ability to balance while in motion or when switching between positions, and both components can be tested to evaluate the risk of falling, through observational tests or measurements (see e.g. [START_REF] Melzer | Postural stability in the elderly: a comparison between fallers and non-fallers[END_REF]).

Static Balance (also called Postural control) is generally considered as the easiest component to measure, as it does not require complex protocols and expensive sensors, and has been shown to be correlated to the risk of fall [START_REF] Melzer | Postural stability in the elderly: a comparison between fallers and non-fallers[END_REF]. Postural control is generally measured by recording the displacement of the center of pressure (CoP) -an approximation of the projection of the Center of Mass, see [START_REF] Benda | Biomechanical relationship between center of gravity and center of pressure during standing[END_REF] -over time using a force platform [START_REF] Piirtola | Force platform measurements as predictors of falls among older people-a review[END_REF]. The resulting signal is usually called a statokinesigram (in this work, a two dimensional signal), that has been shown to contain significant information on postural control [START_REF] Baratto | A new look at posturographic analysis in the clinical context: sway-density versus other parameterization techniques[END_REF] for a wide range of acquisition protocols, such as the Romberg's test [START_REF] Khasnis | Romberg's test[END_REF], [START_REF] Audiffren | A non linear scoring approach for evaluating balance: classification of elderly as fallers and non-fallers[END_REF].

Related Works. The problem of quantifying postural control is hardly new, and numerous previous works have analyzed statokinesigrams with the aim of distinguishing between specific populations. The most commonly used approaches proceed as follows : a small group of scalar features are derived from the CoP trajectory -commonly used features include average speed, lateral variance and surface spanned by the trajectory (see e.g. [START_REF] Baratto | A new look at posturographic analysis in the clinical context: sway-density versus other parameterization techniques[END_REF]). Then a statistical analysis -generally a parametric statistical test -is applied to determine whether those features follow different distributions depending on the populations studied. For instance, in the seminal work of [START_REF] Baratto | A new look at posturographic analysis in the clinical context: sway-density versus other parameterization techniques[END_REF], it is shown that some features such as sway density peaks follow different distributions in parkinsonian, osteoporotic and a heatlhy populations; in [START_REF] Melzer | Postural stability in the elderly: a comparison between fallers and non-fallers[END_REF] others descriptors such as length of the trajectory were shown to be statistically different between fallers and non fallers populations. While these works highlight the effect of different medical conditions on balance, they do not easily translate into reliable postural control quantification method, particularly regarding the risk of falling.

More recently, Machine Learning algorithms (and in particular, scoring algorithms) have been applied to attempt to quantify static balance (see e.g. [START_REF] Audiffren | A non linear scoring approach for evaluating balance: classification of elderly as fallers and non-fallers[END_REF], [START_REF] Liu | A fall detection system using knearest neighbor classifier[END_REF], [START_REF] Bargiotas | On the importance of local dynamics in statokinesigram: A multivariate approach for postural control evaluation in elderly[END_REF]). They show that multivariate approaches (i.e. models using multiple features at the same time, generally through non linear combination) produce significantly better results, and classify fallers and non fallers populations with state of the art performances, with AUC ranging from 0.75 to 0.8. However, these methods relies on algorithms whose learned models are difficult to interpret from a medical point of view, such as Ranking Forest [START_REF] Clémençon | Ranking forests[END_REF] or local aggregation of Gaussian Mixtures. This is the main drawback of these methods, as in medical applications it is important to understand precisely the results in order to A) interpret them to derive new insights regarding the mechanisms of postural control and B) avoid misinterpretation or "mislearning", as errors can have significant consequences [START_REF] Ching | Opportunities and obstacles for deep learning in biology and medicine[END_REF].

Contributions. In this work we introduce a new approach to evaluate postural control with learning algorithms, that achieves strong performances while producing models that are easy to understand and interpret. This approach, detailed in Section 3, relies on • A novel model for postural control and statokinesigrams, from which we derive multiple features,

• A bagging-based algorithm that combines these features while retaining complete interpretability,

• A model-space regularization, based on postural control heuristics, that constrains the search of the optimal classifiers in a restricted manifold and significantly improves the performance of the learning process.

The dataset we used is described in Section 2, and our approach is experimentally studied extensively in Section 4, where we show that this model outperforms its classical counterparts.

Dataset Acquisition

We start by introducing the protocol and measurement device that were used for this study.

Measurements. Data acquisition was performed using a Wii Balance Board (Nintendo, Kyoto, Japan) -a low cost force platform that has gained increasing popularity in Posturography in the recent years [START_REF] Clark | Validity and reliability of the nintendo wii balance board for assessment of standing balance[END_REF]. For each acquisition, we followed the Romberg's test protocol described in [START_REF] Audiffren | A non linear scoring approach for evaluating balance: classification of elderly as fallers and non-fallers[END_REF]: the participant was asked to stand still on the platform 25 seconds eyes opened then 25 seconds eyes closed, resulting in two statokinesigrams : one open-eyes and one closedeyes. The resulting two dimensional signals (containing the medio-lateral and the antero-posterior coordinates of the CoP over time) were then resampled at 25Hz, using the SWARII algorithm [START_REF] Audiffren | Preprocessing the nintendo wii board signal to derive more accurate descriptors of statokinesigrams[END_REF], resulting in signals of length 625 (one example is shown in Figure 1).

Studied populations.

In order to test the relevance of our models, we collected data from three different populations, which were ranked according to prior belief about the qualitative nature of their static balance.

• Group 1 : 67 signals recorded from young healthy individuals aged between 20 and 30. This group was assumed to have the best balance, as postural control tends to worsen with age [START_REF] Tinetti | Preventing falls in elderly persons[END_REF].

• Group 2 : 60 signals recorded from healthy elderly individuals, aged between 74 and 86, which have no medical impairments in balance control and did not report any fall.

• Group 3 : 56 signals recorded from elderly individuals, aged between 74 and 86, which have declared falls in the past 6 months. This group was assumed to have the worst balance due to the report of falls [START_REF] Fuller | Falls in the elderly[END_REF].

It is worth noting that the groups are ordered according to their supposed quality of balance, and that in particular our choice of elderly individual in groups 2 and 3 was in accordance with previous works such as [START_REF] Audiffren | A non linear scoring approach for evaluating balance: classification of elderly as fallers and non-fallers[END_REF], [START_REF] Bargiotas | On the importance of local dynamics in statokinesigram: A multivariate approach for postural control evaluation in elderly[END_REF], thus allowing to compare our performances to previous works. This data acquisition was part of the clinical study registered at ANSM (ID RCB 2014-A00222-45).

Model and Algorithm

As stated in the introduction, in order to evaluate postural control from statokinesigrams we use a novel approach that combines A) a new family of descriptors, that can be easily interpreted while providing an efficient ten dimensional embedding of the statokinesigrams, B) a family of white box learning algorithms based on weak one dimensional classifiers and C) a model-space regularization derived from medical heuristics that greatly improved the performances of the weak classifiers. In this section, we start by briefly introducing the model we chose to derive features from statokinesigrams (Section 3.1), and in the second part we present the families of learning algorithm and we discuss the choice and advantages of a model-based regularization approach (Section 3.2).

Feature Derivation

3.1.1. Posturographic Target. Most of the features we derived from statokinesigrams are inspired by the concept of posturographic target (PT) [START_REF] Baratto | A new look at posturographic analysis in the clinical context: sway-density versus other parameterization techniques[END_REF]. According to this hypothesis, the movement of the CoP during static balance consists of alternating phases of 1) oscillations around a PT and 2) rapid movement to the next PT.

As the nature of the phases and the positions of the PT are unknown, this model is difficult to use in practice. In this work, we consider a different approach with an unique hypothetical posturographic target (HPT) that continuously moves over time -a first approximation to the original model succession of discrete PTs. To compute the position of the HPT, we assume that In our experiments, the average was computed using geometric weights, with more weight being given to the furthest time in the past. As a result, the HPT displacement is not too sensitive to oscillations or momentary variations of the CoP. Although this model is simple, it yielded good performances in our experiments and naturally led to the definition of the following interpretable features.

HPT-derived Features.

From the aforementioned model we derived three groups of features:

• Instability of HPT. These first three features quantify the instability of the HPT, by combining different metrics : the speed of the CoP, the speed of the HPT and the distance between HPT and CoP.

• Excursions. This second group of feature is based on the notion of excursion. We define an excursion as an extended period of time during which the CoP first moves away from the HPT, and then returns close to it. From this concept, we derive the following features -Remoteness duration: time during which CoP is moving away from the HPT. -Return mean time: mean time from which the CoP starts moving closer to the HPT. -Return mean distance: mean distance from which the CoP starts moving closer to the HPT.

• Dispersion Metrics. This third group includes features that are commonly used in Posturography (see e.g. [START_REF] Doyle | Generalizability of center of pressure measures of quiet standing[END_REF]), namely -Ellipse Confidence Area: area of the 95% confidence ellipse containing the entire CoP trajectory -Horizontal dispersion: variance of the horizontal coordinate of the CoP.

In total, we derived 16 features (8 for both the open-eyes and closed-eyes statokinesigrams) from this model, which are used as the signal embedding for the rest of this work.

Scoring Algorithm and Model-space regularization

A key objective of our work was to maximize the interpretability of the learned model while providing accurate evaluation of postural control. This requirement prohibited the use of powerful black box algorithms such as deep learning networks and favored explainable methods such as linear classifiers. Following this idea, we considered the family of learning algorithms based around the combination of easily interpretable weak classifiers (in our case, one dimensional linear support vector machines). These methods include Decision Trees, Adaptive Boosting [START_REF] Freund | A decision-theoretic generalization of on-line learning and an application to boosting[END_REF] and Bagging [START_REF] Breiman | Bagging predictors[END_REF].

Another important aspect of our approach is the choice of regularization, as its appropriate use is of paramount importance for learning algorithms, in particular when applied to medical data. Indeed, such datasets tend to be small in size, and usually contain multiple significant sources of noise, such as protocol related inconsistencies, intra-subject variability (tiredness, ..), sensor noise, etc. As a result, even simple models, that are considered more resilient to overfitting, tend to perform poorly if not properly regularized [START_REF] Schölkopf | Learning with kernels: support vector machines, regularization, optimization, and beyond[END_REF].

Commonly used regularizations include Tykhonov / Lasso regularization, or more generally norm-derived penalties (L p norms, Schatten norms, ...), as well as different approaches such as decision tree pruning. While these regularizations are theoretically sound and powerful off-theshelf tools, in this work we chose to use heuristic based regularizations, that produced significantly better results in our experiments (see Section 4). This process is explained below.

Weak linear classifiers. Let k ∈ {1, . . . , 16} be the index of one of the features defined in section 3.1. Let ν > 0 be the tolerance (a hyperparameter).

Let I = x k i1 , x k i2 , . . . , x k in 1 (resp. J = x k j1 , x k j2 , . . . , x k jn 2
) be the values for this feature for the high-quality static balance group (resp. low-quality). Each weak classifier consists in a one dimensional linear SVM, whose hyperplane τ ∈ R (also called threshold) is chosen as the solution of the following optimization problem:

τ = arg max t∈R v(t)
where v(t) is the value of a threshold and is defined as

v(t) = max {v left (t), v right (t)} (1)
with

v left (t) =          0 if x k ∈J 1 x k < t > n 1 n 2 ν n 2 n 1 x k ∈I 1 x k < t - x k ∈J 1 x k < t otherwise v right (t) =          0 if x k ∈I 1 x k > t > n 2 n 1 ν x k ∈J 1 x k > t - n 2 n 1 x k ∈I 1 x k > t otherwise
represent respectively the left and right values of the threshold. Note that the tolerance ν restricts the number of missclassified individuals below (resp. above) the threshold in v left (resp. v right ). Since the value v is piece-wise constant, the maximum is not unique, and in this case we chose τ following the principle of margin maximization. The value function v varies from the traditional objective function of linear classifiers as it enforces that the values of the feature from the high-quality postural control population should be smaller than the threshold while low-quality static balance related values should be greater -breaking the symmetry of the threshold. By limiting the number of possible classifiers, this constraint acts as a model-space regularization. This particular limitation also reflects the heuristic that, for the features defined in Section 3.1, individuals with better static balance generally present lower values than their worse counterpart. These one dimensional classifiers are used in the three aforementioned learning algorithms: as weak classifiers for both Adaptive Boosting and Bagging, and as node decision rules of Decision Trees.

Heuristic Based Bagging. We also considered a variation of the classical bagging algorithm by using additional heuristic derived constraints. As the open-eyes and the closed-eyes parts of the protocol are different exercices involving different mechanisms of postural control [START_REF] Khasnis | Romberg's test[END_REF], we assumed that they contained complementary information regarding static balance and we constrained the bagging algorithm selection of feature to choose five open-eyes features and five closed-eyes features. A weak classifier is then trained on each of the selected features.

Heuristic Based Decision Trees. In addition to regular decision tree and random forest algorithms, we also used a decision forest algorithm that embedded known posturography heuristics. Following the same idea as for the Heuristic Based Bagging, the forest was constrained to contain two trees, one trained on open-eyes features and one of closed-eyes features. Additionally, we replaced the tree pruning by the following rule: after each split, one of the child node was constrained to be terminal: the left child if the threshold was selected using its v left , or the right child otherwise.

Experimental Evaluations

In this section we present the experiments we conducted to test algorithms introduced in section 3.

Learning Algorithms Comparison

This experiment aimed to observe how the different learning algorithms considered in its paper perform, particularly compared to their more usual counterpart. We evaluated the five following algorithms : In all the experiments we set ν = 0.1, while RF and Adaboost hyperparameters were chosen by cross validation. Each algorithm was evaluated on each pair of groups. When compared together, group 1 and 2 were labeled 0 and 1 respectively, group 1 and 3 were labeled 0 and 1 and group 2 and 3 were labeled 0 and 1 (in this case 0 stands for better postural control and 1 for worse). For each pair of groups, the data was split into train and test (70%-30%). The AUC of the algorithms on the test case was averaged over 50 random train test split.

Results are reported in Table 1. All the algorithms show better performances when the groups to separate are more distant in term of balance (group 1 and group 3). This is coherent with our thought that the models capture the degradation of balance. Overall the algorithms using the heuristic based regularization appear to outperform their classical counterparts, despite RF and Adaboost hyperparameters being tuned with cross validation. This result highlights the importance of proper regularization strategies, and might also suggest that carefully designed constraints, if accessible, are better than off-the-shelf regularization. Finally, it is interesting to note that the algorithm with the best global performance was HBagging, which is also arguably the easiest to interprete as it just compute the average values of the weak classifiers. Figure 2 shows the ROC curves for each pair of group for HBagging.

Comparing Population Scores

This experiment aimed to verify if the quantification of balance computed by HBagging was coherent with the qualitative evaluation of the different groups. To test this hypothesis, we proceeded as follows: group 1 (i.e. individuals with the best postural control) was labeled as class 0 while group 3 (the worst group) was labeled as class 1. Elements of group 1 and group 3 were randomly split into train and test (50%-50%). After training, the model was used to compute the score of each individual in the test set of group 1 and group 3, and in the entire set of group 2. We then compared the mean and standard deviation of each group. This experiment was repeated using 10 different random train test splits.

Results are reported in Table 2. It can be seen that the scores increase with the expected qualitative degradation of balance.

Discussion

Heuristic based models. As illustrated by the results of Table 1, heuristic based learning algorithms and regularizations outperform their regular counterparts while maximizing the resulting model interpretability. This remark encourages further works in this direction.

Score Interpretation. An advantage of Hbagging is the interpretability of the score. In each dimension, the weak classifier provides a class prediction. If the prediction is one, it means that the subject has a statokinesigram classified as abnormal is this dimension. The final score is the average of the scores obtained with the open-eyes classifier and the closed-eyes classifier. Therefore the final score represents the proportion of dimensions with values larger than the one of control population. Higher scores suggest impairments in postural control and therefore increased risk of falling.

Limitation of the model. Our model may not fully and accurately quantify postural control during the Romberg's test. Interpretability is therefore crucial to detect weaknesses in the model in order to A) provide context to the results and B) permit meaningful discussions with medical specialists for continuous adjustements of the features and heuristics. For example in our model the statokinesigram of an individual affected by Parkinson's disease may be given a low score, although it should be classified as abnormal as the individual has an impaired balance. It means that our features and heuristics do not capture yet all the nuances in balance degradation.

Shortcoming in our model include also the lack of local information about the movement of the CoP as all features consist of trajectory-long averages. Therefore, an interesting future research direction would be to adapt our model along the lines developped in [START_REF] Bargiotas | On the importance of local dynamics in statokinesigram: A multivariate approach for postural control evaluation in elderly[END_REF].

It is crucial to note that our model's score alone does not give enough information for a proper diagnosis. In particular it does not take into account eventual particularities of the protocol sequence, for example slight deviation from the theoretical protocol or stress during experiment, or other knowledge about the individual that may change completely the interpretation of the score.

Figure 1 .

 1 Figure 1. Examples of statokinesigrams acquired with the Wii Balance Board, after resampling.

• 1 .

 1 RF -Random Forest• 2. Adaboost -Adaboost with Decision Trees • 3. HBagging -Heuristic based bagging (Section 3) • 4. HDTrees -Heuristic based Decision Trees (Section 3) • 5. HAdaboost -Adaboost with the heuristic based Weak Classifier (Section 3)

Figure 2 .

 2 Figure 2. ROC curves for HBagging

TABLE 1 .

 1 MEAN AND STANDARD DEVIATION OF AUC OF THE DIFFERENT ALGORITHMS FOR EACH PAIR OF GROUPS.

TABLE 2 .

 2 MEAN AND STANDARD DEVIATION OF GROUP SCORES.

Acknowledgments

The authors would like to thank Damien Ricard (COGNAG-G UMR 8257, CNRS, Ecole du Val-de-Grâce, Service de Neurologie, HIA Percy, SSA, France) for its invaluable help in collecting the statokinesigrams. This work was partly supported by the Swiss National Science Foundation under grant number 407540_167320.


Zineb Tabbakh
email: tabbakh.zineb@gmail.com

Mohammed V Lerma-Emi

Ellaia Rachid
email: ellaia@emi.ac.ma

Habbal Abderrahmane
email: habbal@unice.fr

Tabbakh Zineb

New Hybrid Algorithm based on Nonmonotone Spectral Gradient and Simultaneous Perturbation

Keywords: Nonmonotone line search, spectral gradient method, Simultaneous Perturbation Stochastic Approximation

In this paper, we introduce a new hybrid method called Nonmonotone Spectral Gradient and Simultaneous Perturbation (NSGSP). It combines the advantages of nonmonotone spectral gradient (NSG), and simultaneous perturbation (SP) methods. The main idea of our approach is to use the simultaneous perturbation (SP) method in order to get a non expensive estimate of the gradient, and exploit the good properties of the nonmonotone spectral gradient (NSG) method in order to compute an efficient line search. Several numerical experiments are provided. The results indicate that the new method is effective and outperforms most of other popular methods.

is President of the Moroccan Society of Nature-Inspired Computing and Metaheuristics (METAH).

Habbal Abderrahmane is Associate Professor at Nice Sophia-Antipolis University, and INRIA Sophia-Antipolis, France. Research activities of A. Habbal concern analysis and control of systems governed by partial differential equations (PDEs), optimization theory and algorithms and PDE-constrained games. Application fields are related to (nonlinear) mechanics, image processing, data recovering and cell dynamics.

Introduction

We consider minimization problems of the general form minimize f (x) subject to x ∈ S

(1

)
where f is a smooth cost function and S ⊂ R n defines the optimization constraints.

In the vast nonlinear optimization literature, many deterministic and stochastic methods have been developed to solve (1). Among the most popular are the steepest descent algorithms (It is out of the scope of the present paper to provide a review of these methods). From the latter family, spectral gradient methods, originated in [START_REF] Barzilai | Two point step size gradient methods[END_REF] have been applied successfully to find local minimizers of minimization problems for which computational issues arise such as cost and gradient evaluation is expensive, the Hessian is either not available or requires a prohibitive amount of storage. Notice that these methods are essentially steepest descent method, where an exact computation of the optimal step size along the direction of negative gradient is approximated by different effective choices, yielding as many various methods. [START_REF] Barzilai | Two point step size gradient methods[END_REF] step size is computed from a two-point approximation to the secant equation required in quasi-Newton methods. [START_REF] Raydan | The barzilai and borwein gradient method for the large scale unconstrained minimization problem[END_REF] presented the global Barzilai-Borwein (GBB) method based on the technique of nonmonotone line search, which was first introduced in [START_REF] Grippo | A nonmonotone line search technique for newton's method[END_REF]. This technique allows the convergence to a minimum of a broader class of cost functions f , starting from any initial point x 0 , than the usual well behaving case (strongly convex to say it roughly). Its main characteristic is the non-descent property.

Therefore, [START_REF] Birgin | Nonmonotone spectral projected gradient methods on convex sets[END_REF] have improved the work of [START_REF] Raydan | The barzilai and borwein gradient method for the large scale unconstrained minimization problem[END_REF] and introduced the nonmonotone spectral projected gradient (NSPG) method. This approach was developed for the minimization of smooth functions on a closed convex set and has been successfully applied, see e.g. [START_REF] Lakhbab | Hybrid nonmonotone spectral gradient method for the unconstrained minimization problem[END_REF]. On the other hand, simultaneous perturbation stochastic approximation method (SPSA) has been developed by [START_REF] Spall | Adaptive stochastic approximation by the simultaneous perturbation method[END_REF] and used to solve the problem (1). In this method, the gradient is estimated using only two cost evaluations regardless of the dimension n of the optimization space. In simultaneous perturbation stochastic approximation (SPSA), all component directions are perturbed simultaneously using perturbations which are vectors of independent random variables assumed to be symmetric, zero-mean, ±1 valued, and Bernoulli distributed.

Unlike standard finite difference approximations which require 2n cost evaluations to estimate the gradient, [START_REF] Kiefer | Stochastic estimation of the maximum of a regression function[END_REF], SPSA methods are very attractive as soon as the cost evaluation is computationally expensive. Indeed, it is the case for most engineering optimization problems, where cost functions are obtained through the solution of very large scale linear and nonlinear systems, and where closed or implementable forms of the gradient are generally not available.

The main objective of this paper is to propose a method to solve problem (1) combining the nonmonotone spectral gradient (NSG) with the simultaneous perturbation (SP) method. The descent direction is a nearly unbiased estime of the gradient provided by the simultaneous perturbation (SP) method. Along this direction, we use a spectral choice of the step, combined with a nonmonotone line search strategy allowing for the cost function to increase for a limited number of iterations. The advantages of the proposed algorithm are assessed on numerous benchmark functions and the results are compared with many of the existing methods. The obtained results illustrate the efficiency and robustness of our method to solve large scale problems for a broad class of smooth objective functions. The remaining of the paper is organized in the following way. In Section 2, the new hybrid method is described. In Section 3, numerical results are presented for the solution of various benchmark problems. Finally, in Sect.4 the paper is concluded.

NSG, SPSA and hybrid NSGSP methods

We recall in this section the basic features of the Nonmonotone Spectral Gradient (NSG) and Simultaneous Perturbation (SP) methods ; then, we describe with more details our proposed Nonmonotone Spectral Gradient and Simultaneous Perturbation hybrid method (NSGSP), and finally shortly recall the penalization method which we use in (NSGSP) to deal with constrained optimization.

Nonmonotone Spectral Gradient NSG method

The Problem (1) is frequently numerically solved, seeking a local minimizer x * , by using iterative methods according to the following generic scheme

x k+1 = x k -α k g k (x k ) (2)
where g k is the gradient of f at x k , and α k ∈ R + the step length, given initial x 0 ∈ R n .

Instead of an exact and expensive computation of the optimal line search

α k = argmin α f (x k + αg k (x k ))
the so-called spectral gradient methods aim at providing cheap yet effective approximations of the step length. By combining the nonmonotone line search of [START_REF] Grippo | A nonmonotone line search technique for newton's method[END_REF] with the spectral gradient choice of the step length of [START_REF] Barzilai | Two point step size gradient methods[END_REF], [START_REF] Raydan | The barzilai and borwein gradient method for the large scale unconstrained minimization problem[END_REF] proposed the following choice

α k = s T k-1 s k-1 s T k-1 y k-1 (3) where s k-1 = x k -x k-1 and y k-1 = g k -g k-1 .
This line search essentially enforces the following condition:

f (x k+1 ) ≤ max 0≤j≤min{k,M } f (x k-j ) + γg T k (x k+1 -x k ), (4)
where M is a nonnegative integer and γ is a small positive number. Among alternative approaches to compute a line search, let us mention the one based on safeguarded quadratic interpolation proposed in [START_REF] Dai | Modified two-point stepsize gradient method for unconstrained optimization[END_REF]. The following quadratic model is used to approximate the objective function at f (x k + λs k-1 ):

q k (λ) = f k + λg T k s k-1 + 1 2 α k λ 2 ||s k-1 || 2 , ( 5
)
where

α k = s T k-1 y k-1 y T k-1 y k-1
, and the above model satisfies two interpolation conditions ∇q k (0) = g T k s k-1 and q k (-1) = f k-1 , which yields the following choice of stepsize formula

λ k = - s T k-1 s k-1 2(f k-1 -f k ) + 2s T k-1 g k . ( 6
)
After the initial choice of step size, the nonmonotone line search is used to determine the appropriate step length in the appropriate direction, based on the formula (6).

Simultaneous perturbation

In practice, exact values of the gradient of cost functions, notably those considered in engineering applications, are either not available or very costly to approximate by some finite difference scheme, due to the large number of variable and to expensive evaluation of the costs. In this case, nonmonotone spectral gradient (NSG) cannot be applied directly.

To overcome the problem of the possibly prohibitive computational cost in gradient approximation, one may consider the simultaneous perturbation (SP) method due to [START_REF] Spall | Adaptive stochastic approximation by the simultaneous perturbation method[END_REF].

The gradient estimate of f (.) when using simultaneous perturbation (SP) is given by the following centered formula

ĝk (x k ) = f (x k + c k ∆ k ) -f (x k -c k ∆ k ) 2c k ∆ -1 k1 , ∆ -1 k2 , ..., ∆ -1 kn T (7)
where c k is a small positive number, and the components

∆ ki (i = 1, 2, • • • , n) of the perturbation ∆ k = [∆ k1 , ∆ k2 , • • • , ∆ kn ]
T are independent, symmetrically distributed about zero with finite inverse and second moments, and satisfying |∆ ki | < α 0 , and E( |∆ -1 ki ) < α 1 for all indices (k, i) for some finite α 1 0 and α 1 . A simple and popular choice for each component of ∆ k is a Bernoulli ±1 distribution with probability 1 2 . As shown in Spall (1992), simultaneous perturbation (SP) provides an asymptotically unbiased estimate of the exact gradient g k while requiring only a fraction 1 n of the total cost spent in approximating the gradient with a centered finite difference scheme.

We recall hereafter the estimation theorem proved in Spall (1992), as well as a sketched proof (for it is simple enough, despite the striking claim of the theorem).

Theorem 1: Let be ∆ k = [∆ k1 , ∆ k2 , • • • , ∆ kn ]
T that satisfy the above conditions, and suppose that f is smooth enough, at least twice continuously differentiable, and noise-free. Then, ĝk (x k ) is a nearly unbiased estimate of the

g k (x k ) of f E (ĝ k (x k )) = g k (x k ) + O(c 2 k ) almost surely (8)
where g k (x k ) is the true gradient of f at the -deterministic-point x k .

Proof 1: At some fixed iteration k, we take noise-free measurements of f at x k + c k ∆ k , using a standard second order Taylor series expansion,

f (x k ± c k ∆ k ) = f (x k ) ± c k g(x k ) T ∆ k + c 2 k 2 ∆ T k .∇ 2 f (x k )∆ k + O(c 3 k ∆ k 3 ) (9) Then, for i = 1 • • • n, one has E [ĝ ki (x k )] = E f (x k + c k ∆ k ) -f (x k -c k ∆ k ) 2c k ∆ ki (10)
From ( 9) we have

E [ĝ ki (x k )] = E c k g(x k ) T ∆ k + O(c 3 k ∆ k 3 ) --c k g(x k ) T ∆ k + O(c 3 k ∆ k 3 ) 2c k ∆ ki (11)
It follows that:

E [ĝ ki (x k )] = E 2c k g(x k ) T ∆ k + O(c 3 k ∆ k 3 ) 2c k ∆ ki (12) = E     2c k n j=1 g j (x k )∆ kj + O(c 3 k ∆ k 3 ) 2c k ∆ ki     (13) = g i (x k ) + j =i g j (x k )E ∆ kj ∆ ki + O c 2 k E ∆ k 2 . ( 14
)
Hence, under the conditions given above on ∆ k ,

∀j = 1, • • • , n E ∆ 2 kj = 1
, and for all j = i E

∆ kj ∆ ki = 0.
Consequently,

E [ĝ ki (x k )] = g i (x k )(x k ) + O(c 2 k ). ( 15
)
The bias of gradient estimation is of order O(c 2 k ) .

Remark 1: Notice that the proof above assumes that the point x k is deterministic, which no more holds when we shall consider iterates x k updated using the stochastic SP approximation. The proof however still holds in this case, by considering expectations conditional to x k , see Spall (1992).

The proposed NSGSP method

The nonmonotone spectral gradient and simultaneous perturbation (NSGSP) method combines the local search algorithm nonmonotone spectral gradient (NSG) with the descent method where the gradient is estimated by using the simultaneous perturbation. The recursive procedure (2) becomes:

xk+1 = xk -α k ĝk (x k ), ( 16
)
where xk is the current point, ĝk (x k ) is the estimate gradient of f at xk obtained using simultaneous perturbation (SP), and α k is the Barzilai-Borwein (BB) step size attained via nonmonotone line search.

Handling the constraints in NSGSP

To solve constrained optimization problems we have used a classical penalty method, owing to some encouraging properties such as the simplicity of these methods and their ease of implementation. These methods were originally proposed by Courant R. ( 1943) and developed among others by Carroll (1961) and Fiacco and McCormick (1966).

Applying penalty methods, we solve a sequence of unconstrained problems by adding to the objective function a penalty term for violation of constraints. The penalty function is weighted by an increasing sequence of large positive numbers, forcing the unconstrained solution to converge to a feasible one, that is, to satisfy all constraints. There are two types of methods: interior penalty methods, and exterior penalty methods. Since our problem has both inequality and equality constraints, we used exterior penalty functions and combined it with the NSGSP method to handle the constraints in the problem (1).

Assume that the feasible region S of the minimization problem (1) is given by

S = {x : G i (x) ≤ 0, i = 1, 2, • • • , m and H j (x) = 0, j = 1, 2, • • • , p} .
where m is the number of inequality constraints and p the number of equality constraints.

The problem (1) can be replaced by the unconstrained one as follows

min x∈R n f (x) + µP (x), (17)
where µ is a positive parameter and P is a penalty function of the form:

P (x) = m i=1 (max[0, G i (x)]) 2 + p j=1 |H j (x)| 2
P is of course positive and continuous and P (x) = 0 if and only if x ∈ S. However, despite its implementation simplicity, the difficulty of the penalty function approach is to find an appropriate sequence of penalty parameters, that must be chosen large enough to converge to the solution of the constrained problem but not too large, to avoid degeneracy and ill conditioning.

Computational Experiments

The proposed Nonmonotone Spectral Gradient and Simultaneous Perturbation (NSGSP) algorithm was implemented in Matlab. We have followed a Taguchi methodology to tune different parameters specific to the code, and then applied our method to a wide range of benchmark problems. We assessed our algorithm versus a large number of unconstrained and constrained academic problems, and compared our obtained results to numerous popular methods. We thereafter led a comparative study on a worldwide referenced engineering optimization benchmarks : three bar truss, optimal control, pressure vessel design, tension-compression spring design and the speed reducer design. For all the presented numerical experiments, the same stopping criteria were used : either the gradient norm was below some small tolerance threshold (relevant for unconstrained and penalized problems), or the maximum iterations were used without achieving convergence.

Preliminary parameter design and analysis

For the robust tuning of parameters specific to our code, we have used a Taguchi methodology. We have paid particular attention to three parameters : γ, c and M . These parameters are not classical in descent methods, and a sensitivity analysis was then mandatory. After a preliminary analysis of the nonmonotone spectral gradient and simultaneous perturbation (NSGSP) algorithm, we have chosen three levels for each parameter, as described in Table 1. The aim is to estimate the parameters variation effects, in order to set up a robust choice of the latter.

The Taguchi method uses a loss function to measure the deviation between an experimental value and a desired one. See e.g. Peace (1993). This loss function is transformed into a signal-to-noise (S/N) ratio, so that a larger signal-to-noise (S/N) ratio corresponds to a better parameter design, that is, the parameter setting with the highest signal-to-noise (S/N) ratio s the best choice for the experiment. In the present study, the lower-the-better-quality characteristic in signal-to-noise ratio is considered

η = -10 log 1 m m j=1 y 2 j (18)
where y j is the observed data in jth experimentation and m is the number of observations. The L 18 (3 3 ) orthogonal array (refer to Peace (1993) for details) shown in Table 2 was used to select the optimal parameters. After the S/N ratios are measured (Figure 1 and Table 2 ), the next step is to analyze their significance and effects based on analysis of variance (ANOVA). The results are shown in Table 3.

The results of the variance analysis of the parameters design experiment (Table 3) show that the parameter c (B) is the factor that represents the larger effect on the total variation in the function value, with a contribution amount of 37, 91%; the second factor is γ with a contribution of 24,50%. Then follows the effect of M with a contribution of 14,62%.

The Taguchi design experiment resulted in the choice of the following parameters for the nonmonotone spectral gradient and simultaneous perturbation (NSGSP) : γ = 10 -4 , M=10 and c = 0.01.

The simultaneous perturbation (SP) sequence is given by c k = c k τ where τ = 0.101 and k is the iteration counter. The setup for other computational parameters is (see the significance in Algorithm 1) α min = 10 -30 , α max = 10 30 , σ 1 = 0.1, σ 2 = 0.9 and etol = 10 -6 .

Next, we use the following indicators : SD for standard deviation, SR for successful rate. The indicator SR is used to calculate the number of successful runs, that is, when the algorithm generates a desired solution with a required accuracy. We have chosen an accuracy of 10 -4 for all the tests we ran.

Furthermore, to illustrate the statistical significance of the performance of the proposed method, Cohen's d statistical distance is used, for unconstrained optimization examples, to quantify the difference between two means. Each minimization example was run 50 times for twice, and Cohen's d can then be calculated as:

d = Mean 1 -Mean 2 σ pooled where σ pooled = (SD 2 1 + SD 2 2 ) 2 (19)
For our minimization benchmarks, the indicators SD and SR, as well as the Cohen's d statistics are reported in Table 7.

Unconstrained optimization benchmark

We consider a standard, widely used, benchmark list of unconstrained optimization problems. Objective functions, spaces and optima values are listed in Table 4 and in Table 5. These test case functions have different features making them challenging, like non convexity and/or multi-modality. They have known local and global minima.

The NSGSP algorithm introduced in the present paper was ran on these test cases, and the obtained results were compared to numerous other well proved algorithms listed in Table 6. The computational results and statistical metrics are presented in Table 7.

We have plotted in Figure 2 (Eason and Fenton function) and in Figure 3 (Aluffi-Pentini function) the evolution history of the objective functions with respect to the optimization iterations. These two profiles illustrate very clearly the non-descent property of the NSGSP algorithm (inherited from the NSG management of the iterates).

Constrained optimization benchmark

We have considered five constrained optimization test cases as detailed below. In Table 8, we have reported the results obtained with the NSGSP method. These results are compared to the ones available from the literature.

Constrained problem I:

f (x) = (x 1 -2) 2 + (x 2 -1) 2 s.t, g 1 (x) = x 1 -2x 2 + 1 = 0, g 2 (x) = x 2 1 /4 -x 2 2 + 1 ≥ 0, -10 ≤ x 1 , x 2 ≤ 10. (20)
The global minimum is attained at x * = (0.82288, 0.91144), and f (x * ) = 1.3935.

Constrained problem II:

f (x) = (x 2 1 + x 2 -11) 2 + (x 2 1 + x 2 2 -7) 2 s.t, g 1 (x) = 4.84 -(x 1 -0.05) 2 -(x 2 -2.5) 2 ≥ 0, g 2 (x) = x 2 1 + (x 2 -2.5) 2 -4.84 ≥ 0, 0 ≤ x 1 , x 2 ≤ 6. (21)
The global minimum is attained at x * = (2.246826, 2.381865), and f (x * ) = 13.59085, see [START_REF] Bracken | Selected applications of nonlinear programming[END_REF].

Constrained problem III:

f (x) = 5.3578547x 2 3 + 0.8356891x 1 x 5 + 37.293239x 1 -40792.141, s.t, g 1 (x) = 85.334407 + 0.0056858x 2 x 5 + 0.0006262x 1 x 4 -0.0022053x 3 x 5 -92 ≤ 0, g 2 (x) = -85.334407 -0.0056858x 2 x 5 -0.0006262x 1 x 4 + 0.0022053x 3 x 5 ≤ 0, g 3 (x) = 80.51249 + 0.0071317x 2 x 5 + 0.0029955x 1 x 2 + 0.0021813x 2
3 -110 ≤ 0, g 4 (x) = -80.51249 -0.0071317x 2 x 5 -0.0029955x 1 x 2 -0.0021813x 2 3 + 90 ≤ 0, g 5 (x) = 9.300961 + 0.0047026x 3 x 5 + 0.0012547x 1 x 3 + 0.0019085x 3 x 4 -25 ≤ 0, g 6 (x) = -9.300961 -0.0047026x 3 x 5 -0.0012547x 1 x 3 -0.0019085x 3 x 4 + 20 ≤ 0, 78 ≤ x 1 ≤ 102, 33 ≤ x 2 ≤ 45 and 27 ≤ x i ≤ 45(i = 3; 4; 5)..

(22) 3.3.4 Constrained problem IV f (x) = (x 1 -10) 2 + 5(x 2 -12) 2 + x 4 3 + 3(x 4 -11) 2 + 10x 6 5 + 7x 2 6 + x 4 7 -4x 6 x 7 -10x 6 -8x 7 , s.t, g 1 (x) = 127 -2x 2 1 -3x 4 2 -x 3 -4x 2 4 -5x 5 ≥ 0, g 2 (x) = 282 -7x 1 -3x 2 -10x 2 3 -x 4 + x 5 ≥ 0, g 3 (x) = 196 -23x 1 -x 2 2 -6x 2 6 + 8x 7 ≥ 0, g 4 (x) = -4x 2 1 -x 2 2 + 3x 1 x 2 -2x 2 3 -5x 6 + 11x 7 ≥ 0, -10 ≤ x i ≤ 10 (i = 1 • • • 7). ( 23
)
Optimal solution is x * = (2.330499, 1.951372, -0.4775414, 4.365726, -0.6244870, 1.038131, 1.594227), and f (x * ) = 680.6300573.

Constrained problem

V f (x) = x 2 1 + x 2 2 + x 1 x 2 -14x 1 -16x 2 + (x 3 -10) 2 + 4(x 4 -5) 2 + (x 5 -3) 2 + 2(x 6 -1) 2 + 5x 2 7 +7(x 8 -11) 2 + 2(x 9 -10) 2 + (x 10 -7) 2 + 45, s.t, g 1 (x) = 4x 1 + 5x 2 -3x 7 + 9x 8 -105 ≥ 0, g 2 (x) = 10x 1 -8x 1 -17x 7 + 2x 8 ≥ 0, g 3 (x) = -8x 1 + 2x 2 + 5x 9 -2x 10 -12 ≥ 0, g 4 (x) = 3(x 1 -2) 2 + 4(x 2 -3) 2 + 2x 2 3 -7x 4 -120 ≥ 0, g 5 (x) = 5x 2 1 + 8x 2 + (x 3 -6) 2 -2x 4 -40 ≥ 0, g 6 (x) = x 2 1 + 2(x 2 -2) 2 -2x 1 x 2 + 14x 5 -6x 6 ≥ 0, g 7 (x) = 0.5(x 1 -8) 2 + 2(x 2 -4) 2 + 3x 2 5 -x 6 -30 ≥ 0, g 8 (x) = -3x 1 + 6x 2 + 12(x 9 -8) 2 -7x 10 ≥ 0, -10 ≤ x i ≤ 10 (i = 1 • • • 10). ( 24
)
Optimal solution is x * = (2. 1719, 2.3636, 8.7739, 5.9598, 0.9906, 1.4305, 1.3216, 9.8287, 8.2800, 8.3759), and f (x * ) = 24.3062091.

Optimal Control problem

The computational solution of optimal control problems is an important field in optimization theory. It has been studied by many workers, see e.g. [START_REF] Delprat | Control of a parallel hybrid powertrain: Optimal control[END_REF], [START_REF] Delprat | The gradient projection method for solving an optimal control problem[END_REF] and [START_REF] Mehra | A generalized gradient method for optimal control problems with inequality constraints and singular arcs[END_REF].

We consider the following linear quadratic optimal control

min u J(x 0 , u) = ∞ 0 [x(t) T Qx(t) + u(t) T Ru(t)]dt subject to ẋ = Ax + Bu ; x(0) = x 0 ( 25
)
where Q is a positive semidefinite matrix, R is a positive definite matrix, and x 0 is given.

The problem (25) can be solved by using the algebraic Riccati equation (see, for example Broucke (2012)), a powerful tool in linear control theory. Since the algebraic Riccati equation is a nonlinear algebraic equation, solving it analytically in dimensions greater than 3 is intractable in general.

Assumption 1: (A, B) is stabilizable. By stabilizability, there exists a feedback gain K such that the closed loop system

ẋ = (A -BK)x (26)
is stable. The feedback law u = -Kx is clearly admissible. The solution of 26 is given by

x(t) = e (A-BK)t x 0 (27)
The cost criterion is then given by

J(x 0 , -Kx) = x T 0 ∞ 0 e (A-BK) T t (Q + K T RK)e (A-BK)t dtx 0
which is finite. Hence the optimization problem is meaningful.

An NSGSP illustration

Consider the optimal control of a servomotor. Here u(t) is a scalar control and x ∈ R 2 .

   min u J(x 0 , u) = ∞ 0 [x T (t)Qx(t) + ru 2 (t)]dt subject to dx dt = 0 1 0 -1 x + 0 1 u (28)
where Q = 1 0 0 0 , r > 0.

We discretize problem (28) as follow

   minimize h( n+1 1 y 2 i + r n 1 u 2 i ) subject to y i+1 -y i -dt.z i = 0 z i+1 + (-1 + dt)z i -dt.u i = 0 (29)
where n = 20, and x = (y, z); the integration step to approximate the cost function is h = 10 -3 and the ODE time step is dt = 10 -4 .

We choose r = 0.5, u 0 = -5.1 n and x 0 = [3.25; 4.95]. By using nonmonotone spectral gradient and simultaneous perturbation (NSGSP) method in problem (29), we found K ≈ [1.41; 0.96] and J(x 0 , u * ) = 4.5712. The optimal control is given by u * (t) = -1.41y(t) -0.96z(t). The same result is obtained by using the Matlab command lqr.

Pressure vessel design problem

The pressure vessel design problem is a well known benchmark problem introduced in Kannan and Kramer (1994). The main objective is to minimise the cost f (x) of the material mass, forming and welding. There are four design variables: T s(= x 1 ) is the thickness of the shell, T h(= x 2 ) is the spherical head thickness, R(= x 3 ) is the inner radius, and L(= x 4 ) is the length of the cylindrical section. Length unit is in inches.

In the literature, T s and T h are integer multiples of 0.0625, and x 3 , x 4 are continuous variables, such that 40 ≤ x 3 ≤ 80, and 20 ≤ x 4 ≤ 60.

The formulation of this problem can be written as follows:

                   minimize f (x) = 0.6224x 1 x 3 x 4 + 1.7781x 2 x 3 3 + 3.1611x 2 1 x 4 + 19.84x 2 1 x 3 subject to g 1 (x) = 0.0193x 3 -x 1 ≥ 0 g 2 (x) = 0.00954x 3 -x 2 ≥ 0 g 3 (x) = 750.0 × 1728.0 -πx 2 3 x 4 -4 3 πx 2 3 ≥ 0 g 4 (x) = x 4 -240.0 ≥ 0 g 5 (x) = 1.1 -x 1 ≥ 0 g 6 (x) = 0.6 -x 2 ≥ 0 (30)
The NSGSP method was applied to the pressure vessel design problem and the obtained results were compared to those of the literature (De KS Lee et al, Wu, Chow and Sandgren) as illustrated in Table 11.

Tension-compression spring design problem

The tension/compression spring design problem is described in Arora (1989). The objective is to minimize the weight f (x) of the tension compression spring (as shown in Figure 5). The optimization problem involves four nonlinear inequality constraints, namely the minimum deflection, shear stress, surge frequency, outside diameter, and bounds on the design variables. The design variables are wire diameter d(= x1), mean coil diameter D(= x2), and number of active coils P (= x3). The NSGSP obtained solution is given in Table 13 where it is compared to those obtained by other methods.

                             minimize f (x) = (x 3 + 2)x 2 x 2 1 , Subject to g 1 (x) = 1 - x 2 2 x 3 71785x 4 1 ≥ 0, g 2 (x) = 4x 2 2 -x 1 x 2 12566(x 2 x 3 1 -x 4 1 ) + 1 5108x 2 1 -1 ≥ 0, g 3 (x) = 1 - 140.45x 1 x 2 2 x 3 ≥ 0, g 4 (x) = x 2 + x 1 1.5 -1 ≥ 0, 0.05 ≤ x 1 ≤ 2.0, 0.25 ≤ x 2 ≤ 1.3, 2 ≤ x 3 ≤ 15 (31)

Speed reducer design problem

The speed reducer design problem is shown in Figure 6. The objective to minimize is the weight of the speed reducer subject to constraints on on bending stress of the gear teeth, surface stress, transverse deflections of the shafts, stresses in the shaft, and bounds on the design variables. The decision variables are the face width b(= x 1 ), module of teeth m(= x 2 ), number of teeth on pinion z(= x3), length of the first shaft between bearings l 1 (= x 4 ), length of the second shaft between bearings l 2 (= x 5 ), diameter of the first shaft d 1 (= x 6 ), and diameter of the second shaft d 2 (= x 7 ).

The problem statement is the following :

minimize f (x) = 0.7854x 1 x 2 2 (3.3333x 2 3 + 14.9334x 3 -43.0934) -1.508x 1 (x 2 6 + x 2 7 ) +7.4777(x 3 6 + x 3 7 ) Subject to g 1 (x) = 27 x 1 x 2 2 x 3 -1 ≥ 0, g 2 (x) = 397.5 x 1 x 2 2 x 2 3 -1 ≥ 0, g 3 (x) = 1.93x 3 4 x 2 x 3 x 4 6 -1 g 4 (x) = 1.93x 3 5 x 2 x 3 x 4 7 -1 g 5 (x) = (( 745x4 x2x3 ) 2 + 16.9 × 10 6 ) 1/2 110.0x 3 6 -1 ≥ 0, g 6 (x) = (( 745x4 x2x3 ) + 157.5 × 10 6 ) 1/2 85.0x 3 7 -1 g 7 (x) = x 2 x 3 40 -1 g 8 (x) = 5x 2 x 1 -1 g 9 (x) = x 1 12x 2 -1 g 10 (x) = 1.5x 6 + 1.9 x 5 -1 g 11 (x) = 1.1x 7 + 1.9 x 5 -1 2.6 ≤ x 1 ≤ 3.6, 0.7 ≤ x 2 ≤ 0.8, 17 ≤ x 3 ≤ 28, 7.3 ≤ x 4 ≤ 8.3, 7.8 ≤ x 5 ≤ 8.3, 2.6 ≤ x 6 ≤ 3.9, 5.0 ≤ x 7 ≤ 5.5 (32)
The results obtained by NSGSP are presented in Table 15, where they are compared to the results available in the literature.

Three bar truss

The three bar truss problem was introduced by Nowcki(1973). It is often used as a benchmark problem in structural optimization. The objective is to minimize the volume of a statically loaded 3-bar truss, subject to stress constraints on each member. The cross sectional areas A i are within the range 0.1 < A i < 1.

The mathematical programming formulation of this problem can be stated as follows:

                       minimize f (x) = (2 √ 2A 1 + A 2 ) × l subject to g 1 (x) = √ 2A 1 + A 2 √ 2A 2 1 + 2A 1 A 2 P -σ ≤ 0 g 2 (x) = A 2 √ 2A 2 1 + 2A 1 A 2 P -σ ≤ 0 g 3 (x) = 1 A 1 + √ 2A 2 P -σ ≤ 0 0 ≤ A 1 ≤ 1, 0 ≤ A 2 ≤ 1 ( 33
)
where l = 100cm, P = 2KN/cm 2 and σ = 2KN/cm 2 . The NSGSP results are presented in Table 16, and compared to other methods in Table 17.

Conclusion

The present paper introduces a new minimization method called the nonmonotone spectral gradient and simultaneous perturbation (NSGSP). This method is applicable when the gradient of the objective function is not available or is expensive to compute. The gradient is replaced by a simultaneous perturbation stochastic approximation (SPSA). It then applies the Barzilai-Borwein step size and uses the nonmonotone line search techniques, which greatly improves the efficiency of the method. with the non monotone line search, the algorithm was able to escape the local minimum to better one. In many test cases, the new hybrid method exhibited some nice global optimization abilities.

The proposed method was applied to well-known benchmark functions and different classical problems, including twenty unconstrained, five constrained minimization problems, and structural optimization problems. All obtained results contribute to demonstrate the effectiveness of our method. The results were compared to those of published works found in literature. In all examples, nonmonotone spectral gradient and simultaneous perturbation (NSGSP) gave very satisfactory results, which indicates that the nonmonotone spectral gradient and simultaneous perturbation (NSGSP) is competitive, and can be used as an optimization technique for solving engineering optimization problems in challenging (non convex, multimodal) problems. A next step is to adapt the NSGSP approach in order to tackle multiobjective optimization problems, the aim being to capture the Pareto Frontier for large scale problems (those known as many-objective problems). Spall, J. C. (1992). Multivariate stochastic approximation using a simultaneous perturbation gradient approximation.IEEE transactions on automatic control, 37(3), 332-341.

Courant, R. (1943). Variational methods for the solution of problems of equilibrium and vibrations. Bulletin of the American mathematical Society, 49(1), 1-23. Peace, G. S. (1993). Taguchi methods: a hands-on approach. Addison Wesley Publishing Company.

Algorithm 1 Description of the NSG-SP method

Require: x0 ∈ R n , an integer M ≥ 1, a small parameter αmin > 0, a large parameter αmax, a sufficiently decreasing parameter γ ∈ (0, 1) and safeguarding parameters 0 < σ1 < σ2 < 1.

Initially, a, c, and γ, α0 ∈ [αmin, αmax] is arbitrary, a small number , a constant parameter β, and the initial and maximal value of the penalty parameter, respectively r; rmax.

Set x k = x0, f = Fr(x k ), r = r0, set k = 1 1: Apply SP procedure to obtain estimate gradient ĝk 2: Set α k = 1/ ĝ(x0) ∞. 3: while ĝk = 0 do 4: Compute d k = -α k ĝk . 5: Set x = x k + λ k d k . 6: if Pr(x) ≤ and f (x) < f (x k ) then 7:
if Fr(x) < Fr(x k ) then 8:

x k = x 9:
else 10: r = min rmax, βr ; 11:

end if; 12: end if; 13: if f (x) ≤ max 0≤j≤min{k,M } f (x k-j ) + γλg T k d k (34) then 14: define λ k = λ, x k+1 = x, s k = x k+1 -x k , y k = ĝk+1 -ĝk
f 1 = 1 + N i=1 x 2 i 4000 - N i=1 (cos(x i / √ i)) x ∈ R N 0
Rastrigin's function

f 2 = 10N + N i=1 [x 2 i -10 cos(2πx i )]
x ∈ [-5, 12; 5, 12] N 0

Alluffi-Pentiny's function

f 3 = 1 4 x 4 1 -1 2 x 2 1 + 1 10 x 1 + 1 2 x 2 2
x ∈ [-10; 10] 2 -0.352386

Bohachevsky's function

f 4 = x 2 1 + 2x 2 2 -3 10 cos(3πx 1 ) -4 10 cos(4πx 2 ) + 7 10
x ∈ [-100; 100] 2 0 Branin function

f 5 = (x 2 -5.1 4π x 2 1 -5 π x 1 -6) + 10(1 -1 8π ) cos(x 1 ) + 10 x 1 ∈ [-5; 10], x 2 ∈ [0; 15] 0.397887
Camel's function

f 6 = 4x 2 1 -2.1x 4 1 + 1 3 x 5 1 + x 1 x 2 -4x 2 2 + 4x 4 2 x ∈ [-5; 5] 2 -1.0316
Shubert function

f 9 = - 2 i=1 5 j=1 jsin(x i (j + 1)) + 1 x ∈ [-10; 10] 2 -24.06249
Pen Holder Function

f 10 = -exp{-| cos(x 1 ) cos(x 2 )e |1-[(x 2 1 +x 2 2 ) 0.5 /π]| | -1 } x ∈ [-11; 11] 2 -0.963540 f 11 = n i=1 (x 4 i -16x 2 i + 5x i ) x ∈ [-4; 4] n -78.33236
Cosine Mixture Function

f 12 = n i=1 x 2 i -0.1 n i=1 cos(5πx i ) x ∈ [-1; 1] n -0.1n
Rosenbrock Function

f 13 = x 2 1 + 2x 2 2 -0.3 cos(3πx 1 ) -0.4 cos(4πx 2 ) + 0.7 x ∈ R 2 0
Bird Function

f 16 = sin(x 1 ) exp[(1 -cos(x 2 )) 2 ] + cos(x 2 ) exp[(1 -sin(x 1 )) 2 ] + (x 1 -x 2 ) 2 x ∈ [-2π, 2π] 2 -106.764537
Sphere Function

f 17 = n i=1
x 2 i

x ∈ [-5.12, 5.12] n 0 Zakharov Function

f 18 = n i=1 x 2 i + n i=1 0.5 i x i 2 + n i=1 0.5 i x i 4
x ∈ [-100, 100] n 0

Six Hump Camel Back Function

f 19 = (4 -2.1x 2 1 + x 4 1 3 )x 2 1 + x 1 x 2 + (-4 + 4x 2 2 )x 2 2 x 1 ∈ [-3, 3], x 2 ∈ [-2, 2] -1.031600
Exponential function

f 20 = exp(0.5 n i=1 x 2 i ) x ∈ [-1, 1] n 1.00000
Eason and Fenton Function

f 21 = 1 10 12 + x 2 1 + 1 + x 2 2 x 2 1 + x 2 1 x 2 2 + 100 (x 1 x 2 ) 4
x ∈ [0, 10] 2 1.74

f 7 =      0.1 0.2 0.2 0.4 0.4 0.6 0.3          
x ∈ [0; 10] 4 -10.342378

Hartman(H -3, 4) Function 1719, 2.3636, 8.7739, 5.9598, 0.9906, 1.4305, 1.3216, 9.8287, 8.28009, 8.3759 24.30620 Present study 2.1975, 2.3053, 8.7582, 5.0963, 1.0535, 1.6020, 1.3681, 9.8739, 8.2892, 8.204324.314 de KS Lee et al. 2.1552, 2.4076, 8.7780, 5.1020, 0.9676, 1.3576, 1.2877, 9.8004, 8.1878

f 14 = - 4 i=1 c i exp[- 3 j=1 a ij (x j -p ij ) 2 ]

Fiacco

  , A. V., & McCormick, G. P. (1966). Extensions of SUMT for nonlinear programming: equality constraints and extrapolation. Management Science, 12(11), 816-828. Carroll, C. W. (1961). The created response surface technique for optimizing nonlinear, restrained systems. Operations research, 9(2), 169-184.

Figure 1 :

 1 Figure 1: Main Effects Plot for SN ratios for the (M, c, γ) NSGSP parameters Variance source

Figure 2 :Figure 3 :

 23 Figure 2: Eason and Fenton gear inertia problem. Minimum is f (1.74, 2.02) = 1.74Figure 3: Aluffi-Pentini's function. Minimum is f (-1.04, 0) = -0.35

Figure 4 :

 4 Figure 4: Shematic of the pressure vessel design problem

Figure

  Figure 5: Tension-compression spring design

Figure 7 :

 7 Figure 7: Three bar truss

  and go to Step 18.

	Parameters/Level	1	2	3
	M		5	10	20
	c		0.0001	0.01	1
	γ		0.0001	0.01	1
	Table 1 NSGSP parameters and their levels
	M	c		γ	value of f (x)	S/N ratio
	5	0,0001	0,0001	0,3979	8,00452
	5	0,0100	0,0100	0,3979	8,00452
	5	1,0000	1,0000	0,6016	4,41384
	10	0,0001	0,0100	0,3979	8,00452
	10	0,0100	1,0000	0,4386	7,15863
	10	1,0000	0,0001	0,4679	6,59694
	20	0,0001	1,0000	0,3979	8,00452
	20	0,0100	0,0001	0,3979	8,00452
	20	1,0000	0,0100	0,4031	7,89174
	15: else
	16:	λnew ∈ [σ1, σ2λ]			(35)
		set λ = λnew and go to Step 5.
	17: end if
	18: compute b k = s T k y k .
	19: if b k ≤ 0 then
	20:	set α k+1 = αmax
	21: else
	22:	compute a k = s T k s k and α k+1 = min{αmax, max{αmin, a k /b k }}
	23: end if
	24: end while

Table 2

 2 Tuning NSGSP parameters : S/N ratios values of experiments.

Table 3

 3 Results of ANOVA : analysis of variance on the tuned parameters.

Table 4

 4 Test functions

	Test/Benchmark functions

Table 5

 5 Test functions

	Test/Benchmark functions

Table 6

 6 Optimization methods used for performance analysis

	Method	Name (Reference)	Used test function
	PGSL	Probabilistic global search lausanneRaphael and Smith (2003)	f 1
	DE-PSO Differential Evolution-Particle Swarm OptimizationZhang et al. (2009)	f 18
	NHGA	Niche Hybrid Genetic AlgorithmWei and Zhao (2005)	f 18
	ESGAT	simple genetic algorithmRaphael and Smith (2003)	f 1
	Genitor	steady state genetic algorithmRaphael and Smith (2003)	f 1
	pPSA	perturbed Particle Swarm AlgorithmXinchao (2010)	f 17
	GAT	The Genetic Algorithm Toolbox in MATLABYuan et al. (2008)	f 2
	HGA	Hybrid Genetic AlgorithmYuan et al. (2008)	f 2
	MGA	Modified Genetic AlgorithmTsoulos (2008)	f 3 ,f 4 ,f 5 ,f 6 ,f 7 ,f 8 ,f 9
	MACO	The modified ant colony optimizationToksary (2006)	f 13
	SZGA	Successive zooming genetic algorithmKwon et al. (2003)	f 13
	LX-PM	Laplace crossover with Power MutationDeep and Thakur (2007)	f 17
	CTSS	Continuous Tabu Simplex Search	f 18
	GA-PSO Genetic algorithms and particle swarm optimizationKao and Zahara (2008) f 14 ,f 15

Table 7

 7 Computational results of test functions

	Problem no.				NSGSP method				Other Methods
		Dimension	The Best	The Worst	Mean	SD	Median	SR	Cohen d	Method	Obtained result
	1	10	0	2.2204e-15 6.8834e-17 9.7074e-32	0	100% 0.14584	ESGAT	0.0515
										Genitor	0.00496
										PGSL	0.0007
		20	0	4.4409e-16 9.1038e-17 1.3652e-32	0	100% 0.057556	ESGAT	0.0622
										Genitor	0.0240
										PGSL	0.0002
		50	1.1102e-16 1.8874e-15 5.9286e-16 1.5314e-31 4.4409e-16 100%	0	ESGAT	0.0990
										Genitor	0.0170
										PGSL	0.0
	2	10	0	2.8422e-14 5.6843e-16 1.5833e-29	0	100%	0	HGA	1.29e-12
										GAT	1.37
		20	0	2.8422e-14 3.4106e-15 8.5303e-29	0	100% 0.29515	HGA	3.5e-12
										GAT	20.65
		50	0	5.6843e-13 1.2164e-13 1.5123e-26 1.1369e-13 100% -0.10723	HGA	9.09e-12
	3	2	-0.35239	-0.35239	-0.35239	1.4828e-31	-0.35239	100%	0	MGA	-0.352386
	4	2	0	2.2204e-16 8.8818e-18 1.8933e-33	0	100% 0.28577	MGA	0
	5	2	0.39789	0.39789	0.39789	3.469e-21	0.39789	100%	0.2092	MGA	0.397887
	6	2	-1.0316	-1.0316	-1.0316	3.7964e-31	-1.0316	100%	0	MGA	-1.0316
	7	4	-10.1532	-10.1532	-10.1532	3.96e-24	-10.1532	100%	0.1009	MGA	-10.107749
	8	4	-10.4029	-10.4029	-10.4029	2.743e-24	-10.4029	100%	0.0118	MGA	-10.342378
	9	2	-186.7309	-186.7309	-186.7309 1.5304e-13	-186.73	100%	0.245	MGA	-24.06249
	10	2	-1	-0.99998	-1	3.6028e-11	-1	100% -0.13288	Karmakar et al. (2009)	-0.963540
	11	2	-78.3323	-78.332	-78.332	4.2005e-28	-78.332	100%	0	Karmakar et al. (2009)	-78.3323
	12	6	-0.6	-0.6	-0.6	4.3289e-31	-0.6	100%	0	Karmakar et al. (2009)	-0.2
	13	2	3.6778e-16 3.3111e-11 9.3907e-13 2.403e-23 4.8213e-16 100% 0.19809	SZGA	2.98e-8
									MACO	0
	14	3	-3.8634	-3.8634	-3.8634	4.5517e-30	-3.8634	100%	0	GA-PSO	-3.86278
	15	6	-3.3224	-3.3224	-3.3224	9.4663e-32	-3.3224	100%	0	GA-PSO	-3.32237
	16	2	-106.764	-106.76	-106.76	6.7855e-28	-106.76			(Karmakar et al. (2009))	-106.764537
	17	30	6.035e-31	2.35e-26	1.08e-27	1.239e-53	5.493e-29 100%	0.1186	pPSA	5.6e-06
										LX-PM	4.75e-11
	18	30	6.9125e-34 1.0671e-28 4.2324e-30 2.788e-58 7.9539e-32 100% -0.17629	CTSS	3e-07
	19	2	-1.0316	-1.0316	-1.0316	4.1908e-31	-1.0316	100% -0.39693 (Karmakar et al. (2009))	-1.031628
	20	10	1	1	1	4.1612e-31	1	100%	0	(Karmakar et al. (2009))	1.00000

Table 8

 8 Computational results of optimal solution and optimal value test problems

	Problem			NSGSP method
	no.
		Dimension Optimal	The Best	The	Mean	SD
			solution		Worst
	1	10	1.3861	1.3861	1.3861	1.1559e-29 1.3861
	2	2	13.592	13.567	13.532	8.32e-04	13.59
	3	5	-30665	-30666	-30666	3.0882e-05 -30666
	4	7	680.62	680.64	680.68	7.4412e-05 680.64
	5	10	24.305	24.314	24.351	1.60e-05	24.307

Table 9

 9 Computational results of constrained test functions

Table 10

 10 Optimal solution of Pressure vessel design problem

	Design	Present	Sandgren	Wu and Chow De KS Lee et
	variable	method			al.
	x 1	1.1332	1.125	1.125	1.125
	x 2	0.6250	0.625	0.625	0.625
	x 3	58.3267	48.97	58.1978	58.2789
	x 4	43.4925	106.72	44.2930	43.7549
	Cost	7195.8	7980.894	7207.494	7198.433

Table 11

 11 Comparison results of Pressure vessel design problem

Table 12

 12 Optimal solution of Tension-compression spring design

	Methods	The best	The worst	Mean	SD
	NSGSP	0.011868	0.011868	0.011868	3.0093e-36
	GWO	0.0126660	0.0122515	0.0121836	1.085e-05
	CSA	0.0126652	0.0126701	0.0127690	1.357e-06
	CPSO	0.0126747	0.0129240	0.0127330	5.20e-04
	HPSO	0.0126652	0.0127190	0.0127072	1.58e-05

Table 13

 13 Comparison results of the Tension-compression spring design problemTable 14 Optimal solution of Speed reducer design problem

	Figure 6: Speed reducer design

Table 15

 15 Comparison results of Speed reducer design problem

Table 16

 16 Optimal solution of Three bar Truss

	Design	Present	DE -	Ray	and	Tsai	Gandomi et
	variable	method	P CGA	Saini			al.
	A 1	0.7874	0.788675	0.795		0.788	0.78867
	A 2	0.4076	0.408248	0.395		0.408	0.40902
	Weight	263.6789	263.89584	264.3		263.68	263.9716

Table 17

 17 Comparison results of Three bar Truss