# On Metaheuristics for Solving the Parameter Estimation Problem in Dynamic Systems: A Comparative Study

Authors: Gisela C. V. Ramadas; Edite M. G. P. Fernandes; António M. V. Ramadas; Ana Maria A. C. Rocha; M. Fernanda P. Costa
Journal: Journal of Optimization
 (2018)

---

## Abstract

This paper presents an experimental study that aims to compare the practical performance of well-known metaheuristics for solving the parameter estimation problem in a dynamic systems context. The metaheuristics produce good quality approximations to the global solution of a finite small-dimensional nonlinear programming problem that emerges from the application of the sequential numerical direct method to the parameter estimation problem. Using statistical hypotheses testing, significant differences in the performance of the metaheuristics, in terms of the average objective function values and average CPU time, are determined. Furthermore, the best obtained solutions are graphically compared in relative terms by means of the performance profiles. The numerical comparisons with other results in the literature show that the tested metaheuristics are effective in achieving good quality solutions with a reduced computational effort.


---

## Body

## 1. Introduction

The present paper addresses the problem of finding a set of parameter values in a dynamic system model that calibrates the model so that it can reproduce the existing experimental data in the best possible way [1]. This is performed by minimizing an objective function that measures the goodness of the fit. Thus, an objective function, which gives the sum of the squared errors between the model predicted state values and the observed values (at certain time instants of a fixed interval), is minimized, subject to a system of ordinary differential equations (ODE). The dynamic model defined by the system of ODE simulates the time varying processes that take place in a certain time interval. This problem is known in the literature as the dynamic model based parameter estimation (DMbPE) process. Solving this problem is crucial in systems biology and medicine, with a great impact on both pharmaceutical and biotechnological industries [2]. The problem is also very common in the chemical engineering area [3] and it has been extensively used to describe physical phenomena [4]. The DMbPE problem may involve nonlinear differential-algebraic equations and a multiplicity of local solutions may exist. Nonconvexity and ill-conditioning issues of the problem may be addressed with efficient global optimization (GO) methods and proper regularization techniques [5]. The problem is even more difficult than that of an algebraic model. The nonlinear dynamic behavior makes the analytical approach rather complicated, if not impossible, for most of the real phenomena. Numerical direct methods are therefore good alternatives to solve the DMbPE problem. Like any parameter estimation problem, the main advantage is that the intended global minimum is known in advance. The most well-known methods, like Levenberg-Marquardt, gradient descent, and the Nelder-Mead method, which have been extensively used in parameter estimation problems, are local optimization methods and cannot guarantee convergence to a global solution. They are able to exploit a specified search region around a given initial value, but the global solution is difficult to find when the initial value is far from the required solution. Algorithms for GO can be roughly divided into two categories: exact algorithms and heuristics [6].Sörensen and Glover [7] have defined metaheuristic as follows: “A metaheuristic is a high-level problem-independent algorithmicframework that provides a set of guidelines or strategies to develop heuristic optimization algorithms.” Further, an implementation of a heuristic optimization algorithm that follows the strategies laid out in the metaheuristic framework is also referred to as a metaheuristic. This term appeared for the first time in a publication by Glover [8]. Metaheuristics are approximate methods or heuristics that are designed to search for good solutions with less computational effort and time than the more classical algorithms. While heuristics are designed to solve a specific problem, metaheuristics are general-purpose algorithms that can be applied to solve almost any optimization problem. The metaheuristics evaluate a set of potential solutions of an optimization problem and perform a series of operations on those solutions in order to find different and hopefully better solutions. Depending on the way these operations are carried out, a metaheuristic can be classified as a local search, constructive, or population-based metaheuristic. For details, the reader is referred to [7].Most metaheuristics use random procedures that invoke artificial intelligence tools and simulate nature behaviors and their performance does not depend on the properties of the problem to be solved. They are alternative methods to find good approximations to optimal solutions of GO problems. In addition, they are derivative-free and easy to implement [9]. Metaheuristics can locate the vicinity of global solutions with relative efficiency although global optimality cannot be guaranteed in a finite number of iterations. On the other hand, exact methods for GO can guarantee global optimality for certain problems but the required computational effort increases often exponentially with the problem size [10, 11].Previous use of metaheuristics, mainly those that use metaphors based on natural evolution and on behavior of animal swarms, to solve some DMbPE problems showed that they are able to provide good quality solutions when real experimental data and noisy artificial data are considered. In [15], a hybrid metaheuristic algorithm that introduces evolutionary operations, namely, mutation and crossover, into the firefly algorithm (FA), is presented. Two variants of the differential evolution (DE) method—a trigonometric and a modified version—have been implemented in [16, 17]. The authors in [2] developed a new procedure based on the scatter search (SS) methodology. A slightly different version of the SS, which uses the fmincon solver from Matlab (Matlab is a registered trademark of MathWorks, Inc.) as an improvement method, is tested in [12] to solve three well-known DMbPE problems, including the chemical isomerization of α-pinene. This latter problem is analyzed and solved by the FA in [13]. Alternatively, exact methods for GO have also been used to solve DMbPE problems. These include the spatial branch and bound (BB) algorithm to compute a global solution of the partially discretized DMbPE problem [18], the α-BB method to obtain a global solution within an ϵ-precision of the completely discretized DMbPE problem [10], and a method based on interval analysis and on the partially discretized DMbPE [19]. In [14], a deterministic outer approximation approach is applied to the reformulation of the DMbPE problem as a finite NLP by applying a complete discretization using orthogonal collocation on finite elements.Since the problem of estimating the parameters of a dynamic model is important, the contribution of this study is concerned with the implementation and practical comparison of five very simple and easy to implement metaheuristics, hybridized with a local intensification phase, when solving the finite nonlinear programming (NLP) problem that arises when a numerical direct method of a sequential type is used to locate a global optimal solution to the DMbPE problem. The unknown parameters are the decision variables of the NLP problem.The selected metaheuristics are very popular and have been used to solve a variety of real-world applications. The selection includes the FA [20], the Harmony Search (HS) algorithm [21], the DE algorithm [22], the Artificial Bee Colony (ABC) algorithm [23], and the Directed Tabu Search (DTS) algorithm [24]. During the local intensification phase, the well-known Hooke-and-Jeeves (HJ) local search [25] is implemented with all the metaheuristics. These local search enhancements lead to much faster convergence and provide good quality solutions. The practical comparison carried out in the present study aims to analyze the performance of each metaheuristic in terms of the quality of the obtained solutions. To test the five metaheuristics, nine DMbPE problems were selected from the literature [4, 10, 13, 19, 26], yielding 12 instances due to different experimental data. The problems are described in the Appendix.The remainder of the paper proceeds as follows. Section2 addresses the DMbPE problem and the sequential numerical direct method, and Section 3 introduces the selected metaheuristics. Section 4 contains the results of all the numerical experiments and the conclusions are summarized in Section 5.


## 2. Sequential Direct Method to Solve the DMbPE Problem

To solve the DMbPE problem, the sequential numerical direct method is used. For completeness, the DMbPE problem reads as follows. Findp∗ such that(1)Jp∗=minpJp≡∑j=1D∑i=1Oyjti-yj,iobs2subject  toy′t=gt,yt,p∀t∈t0,tfyt0=y10,…,yD0pL≤p≤pU,where J:Rn→R represents the objective function that depends on p=(p1,p2,…,pn), a vector of the n decision variables (the parameters that need to be estimated); y=(y1,y2,…,yD) is the vector (with length D) of the state variables that depend on t (the independent variable); and y′ represents the vector with the derivatives of y with respect to t. O represents the number of experimental data available for each yj; yj,iobs, j=1,…,D, i=1,…,O, are the experimentally observed values; y′=g(t,y,p) is the system of first-order ODE; y10,…,yD0 is the vector of initial values to the variables yj, j=1,…,D; and pL and pU represent the vectors with the lower and upper bounds on the parameters, respectively. t0 and tf are the initial and final values, respectively, of the independent variable t.To solve problem (1), there are indirect methods and direct methods [27]. Indirect methods use the first-order necessary conditions from Pontryagin’s minimum principle to reformulate the problem as a two-point boundary value problem. The latter may become difficult to solve especially if the problem contains state variable constraints.In direct methods, the optimization present in (1) is performed directly. This type of method discretizes the problem and applies NLP techniques to the resulting finite-dimensional optimization problem. There are two types of direct methods. In thesequential direct method, the DMbPE problem is transcribed into a small finite-dimensional optimization problem through the discretization of the decision variables only. Hence, the optimization is carried out in the space of the decision variables only. Given a set of values for the decision variables, the system of ODE is accurately integrated (over the entire time interval) using a specific numerical integration formula (with error control mechanisms to enforce the accuracy of the state variable values), so that the objective function value can be evaluated [3]. Thus, the ODE are satisfied at each iteration of the NLP solver. The method is called sequential, since the processes of minimizing the sum of the squared errors and solving the ODE are done in a sequential manner. This method may lead to a slow convergence process since the system of ODE is solved again and again each time an objective function value is required. In thesimultaneous direct method, approximations to the solution of the system of ODE are computed instead. The optimization is carried out in the full space of the discretized decision and state variables, which may lead to very-large-dimensional finite NLP problem with equality constraints (emerging from the discretization of the system of ODE), in particular if a fine grid of points is required to obtain a high level of integration accuracy. Thus, very efficient NLP solvers are required when solving the finite problem.We use the sequential direct method and the system of ODE is numerically integrated by the Matlab functionode15s. The default value for the scalar relative error tolerance is set to 1.0E-08 and all the components of the vector of absolute error tolerance are set to 1.0E-08. The sequentially discretized NLP problem is then solved by a stochastic metaheuristic. This is further elaborated in the next section.


## 3. Stochastic Metaheuristics

To compute global optimal solutions to NLP problems, stochastic or deterministic methods are available. Stochastic GO methods are able to provide a near-optimal solution in a short CPU time, although it may not be globally optimal. In contrast, deterministic GO methods provide an interval within which the global optimal solution falls, although they require very large computational efforts [14]. Stochastic metaheuristics, such as FA, HS, DE, ABC, and DTS, include randomly generated numbers from probability distributions to define a set of solutions on the search region, as well as to move those solutions to hopefully better positions. They ensure convergence to a global solution in probability, while deterministic methods tend to guarantee asymptotic convergence. Global search techniques use exploration and exploitation search procedures that aim to diversify the search so that a global optimal solution is located and to intensify the search so that a good approximation is computed in a promising area of the search space.Most stochastic metaheuristics are classified in terms of source of inspiration as nature-inspired algorithms. The well-known swarm-intelligence-based algorithms belong to a wider class called the bioinspired algorithms, and these are a subclass of nature-inspired algorithms [28]. The most popular are the algorithms based on swarm intelligence, with the FA and the ABC algorithm being two examples. On the other hand, the HS is not a bioinspired algorithm and has its inspiration from music. Although the DE algorithm is not based on any biological behavior, it can be classified as a bioinspired algorithm due to the keyword “evolution,” but it is not swarm-intelligence-based [28]. Finally, the DTS algorithm is not a nature-inspired algorithm, but it has the particularity of using past information of the search to record it in memory. The majority of the metaheuristics define a set of approximate solutions at each iteration, called population, and generate new trial solutions that may be accepted as the solutions for the next iteration, if some improvement is detected relative to the current set of solutions.We use the notationpi=(p1i,p2i,…,pni)T to represent the ith point of a population of NP points, possible solutions of problem (1). The set Ω={p∈Rn:pjL≤pj≤pjU,j=1,…,n} defines the bound constraints of the problem.


### 3.1. The FA

The FA is a bioinspired metaheuristic algorithm that is capable of converging to a global solution of an optimization problem. It is inspired by the flashing behavior of fireflies at night [20, 29]. The main rules used to construct the FA are as follows: (i) all fireflies are unisex, meaning that any firefly can be attracted to any other brighter one; (ii) the brightness of a firefly is determined from the encoded objective function; (iii) attractiveness is directly proportional to brightness but decreases with distance.Table1 lists the most relevant nomenclature used in the basic FA. At the beginning, the positions of the fireflies in the search space are randomly generated. Then, they are evaluated by computing the corresponding objective function value and the points are ranked in ascending order; that is, p1 is the brightest firefly, p2 is the second brightest, and so on. In the classical FA, a firefly i (i=2,…,NP) is moved towards the brighter fireflies j=1,…,i-1 as follows:(2)pi=pi+β0-βmin⁡exp⁡-γpi-pj2+βmin⁡pj-pi+αλi-0.5X,where the second term on the right-hand side of (2) is due to the attraction while the third term is due to randomization with α∈(0,1), β0 is the attraction parameter when the distance is zero, λi is a number uniformly distributed in [0,1], and X is a vector that aims to scale the movement to the set Ω. The parameter γ is crucial to speed the convergence of the algorithm and it can take any value in the set [0,∞). When γ→∞, the attractiveness tends to the minimum constant value βmin. To accelerate convergence and at the same time prevent premature convergence, the sequence of α values is made to slowly decrease along the iterative process. If a component j of the final position of a firefly falls outside the set [pjL,pjU], it is shifted onto the boundaries.Table 1
Nomenclature for the FA.
N
P
Number of fireflies in the population
β
Attractiveness of a firefly
γ
Absorption coefficient/variation of attractiveness
α
Randomization parameter
p
1
The brightest firefly
p
N
P
The less brighter firefly


### 3.2. The HS Algorithm

The HS algorithm was developed to solve GO problems in an analogy with the music improvisation process where music players improvise the pitches of their instruments to obtain better harmony [21, 30]. At each iteration k, the basic HS algorithm provides a set of solution vectors, available in the harmony memory (HM), from which the best and the worst solutions, in terms of objective function values, are identified. Table 2 lists the most relevant nomenclature used in the HS algorithm.Table 2
Nomenclature for the HS algorithm.
HM
Harmony memory
H
M
S
Size of the HM
p
best
The best solution in HM
p
worst
The worst solution in HM
H
M
C
R
Harmony memory considering rate
P
A
R
The pitch adjusting rate
B
W
Distance bandwidthThe HM containsHMS solution vectors that are maintained in memory throughout the iterative process. First, the solutions in the HM are randomly generated in the search space Ω, pi, i=1,…,HMS. Then, they are evaluated and the best harmony, pbest, and the worst, pworst, in terms of objective function value are identified. Hereafter, at each iteration, a new harmony is improvised; that is, a new vector v is generated, using three improvisation operators. When theHM operator is used, the component j of v is chosen from the HM with probability HMCR; otherwise, therandom selection operator is used and the component is randomly generated in Ω:(3)vj=pji,irandom∈1,…,HMS,ifλ<HMCRpjL+λpjU-pjL,otherwise.for j=1,…,n, where λ is a random number uniformly distributed in [0,1]. Finally, thepitch adjustment operator is subsequently applied with a probability 0<PAR(K)<1, which varies with the iteration counter K, to refine only the components j produced by theHM operator, as follows:(4)vj=vj±λBWK,ifλ<PARKvj,otherwise,where 0<BW(K)≤1 is the distance bandwidth that depends also on K [31]. Finally, the components of v are checked against the bounds and projected onto the boundaries if they fall outside. In the final stage, if f(v)<f(pworst), the HM is updated since the new harmony is included in the HM, replacing the worst one.


### 3.3. The DE Algorithm

The DE is a bioinspired population-based algorithm that relies on three strategies—mutation,crossover, andselection—to define the NP solutions/points for the next iteration [22]. The most important nomenclature and parameters are shown in Table 3.Table 3
Nomenclature for the DE algorithm.
N
P
Number of points in the population
F
Amplification parameter
v
∈
R
n
Mutant point
u
∈
R
n
Trial point
C
R
Crossover parameterThe initial population of points,pi∈Rn, i=1,…,NP, is randomly generated in the search space Ω. The most commonly usedmutation strategy defines the mutant point i, vi, as follows:(5)vi=pr1+Fpr2-pr3with uniformly chosen random indices r1, r2, and r3 from the set {1,2,…,NP}, mutually different, and F is a real parameter in [0,2] which controls the amplification of the differential variation, pr2-pr3. The indices r1, r2, and r3 are also chosen to be different from the index i. pr1 is called the base point.The components of the mutant vector are then mixed with components of thep vector to generate the trial point, ui. This strategy is referred to ascrossover and can be described as follows:(6)uji=vji,ifλ≤CRorj=sipji,otherwise.for j=1,…,n, where λ denotes a random number uniformly generated within the interval [0,1] and aims to perform the mixing of the component j of the points, CR∈[0,1], and si, an index randomly selected from {1,…,n}, aims to ensure that ui gets at least one component from vi. Finally, the components of u are projected onto the boundaries of [pL,pU] if they fall outside. Then, aselection strategy compares each trial point ui with pi in terms of objective function values and the best one is selected for the population of the next iteration. The variant described above is referred to as DE/rand/1/bin, where “rand” specifies the vector that is mutated, “1” defines the number of difference vectors, and “bin” means that crossover is based on independent binomial experiments [22]. There are other frequently used DE variants, for instance, the DE/best/1/bin which uses the best point of the population as the base point to define the mutant point and the DE/rand-to-best/1/bin with two differential variations:(7)vi=pi+Fpbest-pi+Fpr1-pr2,where pbest is the best point in the current population [32].


### 3.4. The ABC Algorithm

The ABC algorithm is an optimization algorithm based on the intelligent behavior of honeybee swarms [23]. The colony of artificial bees of size NP includes the employed bees, the onlooker bees, and the scout bees. The first half of the colony consists of employed bees and the second half constitutes the onlookers and scouts. The position of a food source represents a possible solution of the optimization problem and the amount of nectar in the food source gives the quality of that solution. The number of food sources, Nf, is taken to be equal to the number of employed bees. The most important nomenclature and parameters in the ABC algorithm are shown in Table 4.Table 4
Nomenclature for the ABC algorithm.
N
P
Number of bees in the colony
N
f
Number of food sources
p
∈
R
n
Food source/solution
v
∈
R
n
Mutant solution
“limit”
Limit for abandonmentAt the initial stage, a set of food source positions are randomly selected by the bees; that is, the positions are randomly generated in the search spaceΩ, and their nectar amounts are determined in terms of fitness values. During the employed bee phase, new candidate food positions, called mutant solutions, are produced as follows:(8)vji=pji+λpji-pjk,for i=1,…,Nf, where k (different from i) and j are randomly chosen indexes from the sets {1,…,Nf} and {1,…,n}, respectively, and λ is a random number uniformly generated within the interval [-1,1]. If the mutant component vji falls outside [pjL,pjU], it is shifted to the boundaries. Afterwards, the mutant solution vi is compared with pi and a greedy selection is applied to choose the one with better nectar, that is, with better fitness, as described in (9). If the current solution pi is maintained (meaning that it has not been improved), its trial counter is increased.On the other hand, during the onlooker bee phase, the food sources are randomly chosen according to probability values,Pi (associated with the food sources), that depend on their fitness Fi evaluated by(9)Fi=Jpi+1-1,ifJpi≥01+Jpi,otherwise.The mutant solution is produced from the chosen old one as previously shown in (8) and the greedy selection is applied between the current pi and its mutant. Finally, in the scout bee phase, if the trial counter of the solution that has not been improved mostly exceeds a predetermined value, called “limit,” that solution is abandoned and replaced by a position randomly generated in the search space Ω [33].


### 3.5. The DTS Algorithm

The tabu search (TS) algorithm, introduced to continuous optimization in the paper [34], is capable of guiding the search out of local optima and exploring new regions for the global optimum. This is a point-by-point iterative procedure that maintains a list of the most recent movements, denoted by “tabu list” (TL), so that the movements that lead to solutions previously visited are avoided.The DTS method developed in [24] is hybridization of the TS with a direct search method that aims to search the neighborhood of a local minimum. The method comprises three main search procedures: exploration, diversification, and intensification. The main loop consists of the exploration and the diversification search procedures. The exploration search aims to explore the search space Ω and uses direct search methods to be able to stabilize the search, in particular in the vicinity of a local minimum. Cycling is prevented by the TL and by another four TS memory elements: the multiranked tabu list, the tabu region, the semitabu region, and the visited region list (VRL). The VRL works as a diversification tool and is used to direct the search towards regions that have not been visited in the search space. The most important nomenclature and parameters of the DTS algorithm are listed in Table 5.Table 5
Nomenclature for the DTS algorithm.
TL
Tabu list
VRL
Visited region list
H
Ratio of accepting diversification point
E
d
g
=
f
edg
m
i
n
j
(
p
j
U
-
p
j
L
)
Edge length (with0<fedg<1)
d
V
R
=
r
d
V
R
E
d
g
Region radius in VRL (rdVR>0)During the exploration procedure, DTS uses an adaptive pattern search strategy to generate an approximate descent direction (ADD) for the objective functionJ at the current approximation. First, based on an approximation p¯, n pattern directions parallel to the coordinate axes are constructed and n trial points xi, i=1,…,n, are generated, along these directions with a given step length. Thus, the ADD v is obtained by(10)v=∑i=1nwiuiwherewi=Jxi-Jp¯∑j=1nJxj-Jp¯,ui=-xi-p¯xi-p¯.Second, two local trial points xn+1 and xn+2 are then generated along v with two different step lengths aiming to explore the region along v [24]. The diversification procedure aims to randomly generate a new trial point p~∈Ω outside the previously visited regions. The VRL contains the centers of the visited regions and the frequency with which these regions are visited. If the shortest distance of p~ to a region’s center, weighted by its frequency, exceeds the predefined region radius, given by HdVR, the point is accepted; otherwise, a new trial point is generated. When one of the best obtained trial solutions is sufficiently close to a global minimum, or its value has not been changed for a specified number of iterations, the DTS algorithm leaves the exploration and diversification search procedures and invokes an intensification procedure. Here, the HJ local search is used to compute a solution still closer to the global minimum (see the next section).


### 3.6. Hooke-and-Jeeves Local Search

A pattern search method directs the search towards a minimizer using a pattern of specific number points. At leastn+1 points must be provided by the pattern, where n is the number of variables. Based on a current approximation, the HJ method uses a pattern search strategy to define a trial solution that gives an improvement in the objective function value, when compared with the current point [25]. Let the current approximation be denoted by p(ki), where ki represents the iteration counter in this iterative process. Then, a trial point, vi, is generated along a search direction (starting from the current point) with a certain step size Δ(ki)>0 as follows:(11)vi=pki+Δkidi,where di is the search direction chosen from a finite set D of positive spanning directions in Rn. The most used set contains the 2n coordinate directions, defined as the positive and negative unit coordinate vectors D=e1,…,en,-e1,…,-en. The most important property of this set is that at least one of the coordinate directions is a descent direction for the objective J, when p(ki) is not a stationary point of J. When the search fails to generate a trial point that is better than p(ki), the iteration is calledunsuccessful, the step size Δ(ki) is halved so that a refined search can be carried out, and p(ki+1)←p(ki). If Δ(ki) falls below a given stopping tolerance, the algorithm terminates and the current p(ki) is considered an approximate minimizer of J. However, when at the end of each iteration the objective function value has reduced, then the iteration is calledsuccessful, Δ(ki) is not changed, and p(ki+1)←vi.For a fair comparison, and to improve the quality of the produced solutions, we propose the implementation of the HJ intensification phase with all the above-mentioned metaheuristics. Based on the final solution provided by the metaheuristic, the HJ local search algorithm is invoked and allowed to run for5n iterations. The number of function evaluations required by the HJ local search is added to the function evaluations of the metaheuristic, herein denoted by NFEval, to produce the total number of evaluations of the run.


## 3.1. The FA

The FA is a bioinspired metaheuristic algorithm that is capable of converging to a global solution of an optimization problem. It is inspired by the flashing behavior of fireflies at night [20, 29]. The main rules used to construct the FA are as follows: (i) all fireflies are unisex, meaning that any firefly can be attracted to any other brighter one; (ii) the brightness of a firefly is determined from the encoded objective function; (iii) attractiveness is directly proportional to brightness but decreases with distance.Table1 lists the most relevant nomenclature used in the basic FA. At the beginning, the positions of the fireflies in the search space are randomly generated. Then, they are evaluated by computing the corresponding objective function value and the points are ranked in ascending order; that is, p1 is the brightest firefly, p2 is the second brightest, and so on. In the classical FA, a firefly i (i=2,…,NP) is moved towards the brighter fireflies j=1,…,i-1 as follows:(2)pi=pi+β0-βmin⁡exp⁡-γpi-pj2+βmin⁡pj-pi+αλi-0.5X,where the second term on the right-hand side of (2) is due to the attraction while the third term is due to randomization with α∈(0,1), β0 is the attraction parameter when the distance is zero, λi is a number uniformly distributed in [0,1], and X is a vector that aims to scale the movement to the set Ω. The parameter γ is crucial to speed the convergence of the algorithm and it can take any value in the set [0,∞). When γ→∞, the attractiveness tends to the minimum constant value βmin. To accelerate convergence and at the same time prevent premature convergence, the sequence of α values is made to slowly decrease along the iterative process. If a component j of the final position of a firefly falls outside the set [pjL,pjU], it is shifted onto the boundaries.Table 1
Nomenclature for the FA.
N
P
Number of fireflies in the population
β
Attractiveness of a firefly
γ
Absorption coefficient/variation of attractiveness
α
Randomization parameter
p
1
The brightest firefly
p
N
P
The less brighter firefly


## 3.2. The HS Algorithm

The HS algorithm was developed to solve GO problems in an analogy with the music improvisation process where music players improvise the pitches of their instruments to obtain better harmony [21, 30]. At each iteration k, the basic HS algorithm provides a set of solution vectors, available in the harmony memory (HM), from which the best and the worst solutions, in terms of objective function values, are identified. Table 2 lists the most relevant nomenclature used in the HS algorithm.Table 2
Nomenclature for the HS algorithm.
HM
Harmony memory
H
M
S
Size of the HM
p
best
The best solution in HM
p
worst
The worst solution in HM
H
M
C
R
Harmony memory considering rate
P
A
R
The pitch adjusting rate
B
W
Distance bandwidthThe HM containsHMS solution vectors that are maintained in memory throughout the iterative process. First, the solutions in the HM are randomly generated in the search space Ω, pi, i=1,…,HMS. Then, they are evaluated and the best harmony, pbest, and the worst, pworst, in terms of objective function value are identified. Hereafter, at each iteration, a new harmony is improvised; that is, a new vector v is generated, using three improvisation operators. When theHM operator is used, the component j of v is chosen from the HM with probability HMCR; otherwise, therandom selection operator is used and the component is randomly generated in Ω:(3)vj=pji,irandom∈1,…,HMS,ifλ<HMCRpjL+λpjU-pjL,otherwise.for j=1,…,n, where λ is a random number uniformly distributed in [0,1]. Finally, thepitch adjustment operator is subsequently applied with a probability 0<PAR(K)<1, which varies with the iteration counter K, to refine only the components j produced by theHM operator, as follows:(4)vj=vj±λBWK,ifλ<PARKvj,otherwise,where 0<BW(K)≤1 is the distance bandwidth that depends also on K [31]. Finally, the components of v are checked against the bounds and projected onto the boundaries if they fall outside. In the final stage, if f(v)<f(pworst), the HM is updated since the new harmony is included in the HM, replacing the worst one.


## 3.3. The DE Algorithm

The DE is a bioinspired population-based algorithm that relies on three strategies—mutation,crossover, andselection—to define the NP solutions/points for the next iteration [22]. The most important nomenclature and parameters are shown in Table 3.Table 3
Nomenclature for the DE algorithm.
N
P
Number of points in the population
F
Amplification parameter
v
∈
R
n
Mutant point
u
∈
R
n
Trial point
C
R
Crossover parameterThe initial population of points,pi∈Rn, i=1,…,NP, is randomly generated in the search space Ω. The most commonly usedmutation strategy defines the mutant point i, vi, as follows:(5)vi=pr1+Fpr2-pr3with uniformly chosen random indices r1, r2, and r3 from the set {1,2,…,NP}, mutually different, and F is a real parameter in [0,2] which controls the amplification of the differential variation, pr2-pr3. The indices r1, r2, and r3 are also chosen to be different from the index i. pr1 is called the base point.The components of the mutant vector are then mixed with components of thep vector to generate the trial point, ui. This strategy is referred to ascrossover and can be described as follows:(6)uji=vji,ifλ≤CRorj=sipji,otherwise.for j=1,…,n, where λ denotes a random number uniformly generated within the interval [0,1] and aims to perform the mixing of the component j of the points, CR∈[0,1], and si, an index randomly selected from {1,…,n}, aims to ensure that ui gets at least one component from vi. Finally, the components of u are projected onto the boundaries of [pL,pU] if they fall outside. Then, aselection strategy compares each trial point ui with pi in terms of objective function values and the best one is selected for the population of the next iteration. The variant described above is referred to as DE/rand/1/bin, where “rand” specifies the vector that is mutated, “1” defines the number of difference vectors, and “bin” means that crossover is based on independent binomial experiments [22]. There are other frequently used DE variants, for instance, the DE/best/1/bin which uses the best point of the population as the base point to define the mutant point and the DE/rand-to-best/1/bin with two differential variations:(7)vi=pi+Fpbest-pi+Fpr1-pr2,where pbest is the best point in the current population [32].


## 3.4. The ABC Algorithm

The ABC algorithm is an optimization algorithm based on the intelligent behavior of honeybee swarms [23]. The colony of artificial bees of size NP includes the employed bees, the onlooker bees, and the scout bees. The first half of the colony consists of employed bees and the second half constitutes the onlookers and scouts. The position of a food source represents a possible solution of the optimization problem and the amount of nectar in the food source gives the quality of that solution. The number of food sources, Nf, is taken to be equal to the number of employed bees. The most important nomenclature and parameters in the ABC algorithm are shown in Table 4.Table 4
Nomenclature for the ABC algorithm.
N
P
Number of bees in the colony
N
f
Number of food sources
p
∈
R
n
Food source/solution
v
∈
R
n
Mutant solution
“limit”
Limit for abandonmentAt the initial stage, a set of food source positions are randomly selected by the bees; that is, the positions are randomly generated in the search spaceΩ, and their nectar amounts are determined in terms of fitness values. During the employed bee phase, new candidate food positions, called mutant solutions, are produced as follows:(8)vji=pji+λpji-pjk,for i=1,…,Nf, where k (different from i) and j are randomly chosen indexes from the sets {1,…,Nf} and {1,…,n}, respectively, and λ is a random number uniformly generated within the interval [-1,1]. If the mutant component vji falls outside [pjL,pjU], it is shifted to the boundaries. Afterwards, the mutant solution vi is compared with pi and a greedy selection is applied to choose the one with better nectar, that is, with better fitness, as described in (9). If the current solution pi is maintained (meaning that it has not been improved), its trial counter is increased.On the other hand, during the onlooker bee phase, the food sources are randomly chosen according to probability values,Pi (associated with the food sources), that depend on their fitness Fi evaluated by(9)Fi=Jpi+1-1,ifJpi≥01+Jpi,otherwise.The mutant solution is produced from the chosen old one as previously shown in (8) and the greedy selection is applied between the current pi and its mutant. Finally, in the scout bee phase, if the trial counter of the solution that has not been improved mostly exceeds a predetermined value, called “limit,” that solution is abandoned and replaced by a position randomly generated in the search space Ω [33].


## 3.5. The DTS Algorithm

The tabu search (TS) algorithm, introduced to continuous optimization in the paper [34], is capable of guiding the search out of local optima and exploring new regions for the global optimum. This is a point-by-point iterative procedure that maintains a list of the most recent movements, denoted by “tabu list” (TL), so that the movements that lead to solutions previously visited are avoided.The DTS method developed in [24] is hybridization of the TS with a direct search method that aims to search the neighborhood of a local minimum. The method comprises three main search procedures: exploration, diversification, and intensification. The main loop consists of the exploration and the diversification search procedures. The exploration search aims to explore the search space Ω and uses direct search methods to be able to stabilize the search, in particular in the vicinity of a local minimum. Cycling is prevented by the TL and by another four TS memory elements: the multiranked tabu list, the tabu region, the semitabu region, and the visited region list (VRL). The VRL works as a diversification tool and is used to direct the search towards regions that have not been visited in the search space. The most important nomenclature and parameters of the DTS algorithm are listed in Table 5.Table 5
Nomenclature for the DTS algorithm.
TL
Tabu list
VRL
Visited region list
H
Ratio of accepting diversification point
E
d
g
=
f
edg
m
i
n
j
(
p
j
U
-
p
j
L
)
Edge length (with0<fedg<1)
d
V
R
=
r
d
V
R
E
d
g
Region radius in VRL (rdVR>0)During the exploration procedure, DTS uses an adaptive pattern search strategy to generate an approximate descent direction (ADD) for the objective functionJ at the current approximation. First, based on an approximation p¯, n pattern directions parallel to the coordinate axes are constructed and n trial points xi, i=1,…,n, are generated, along these directions with a given step length. Thus, the ADD v is obtained by(10)v=∑i=1nwiuiwherewi=Jxi-Jp¯∑j=1nJxj-Jp¯,ui=-xi-p¯xi-p¯.Second, two local trial points xn+1 and xn+2 are then generated along v with two different step lengths aiming to explore the region along v [24]. The diversification procedure aims to randomly generate a new trial point p~∈Ω outside the previously visited regions. The VRL contains the centers of the visited regions and the frequency with which these regions are visited. If the shortest distance of p~ to a region’s center, weighted by its frequency, exceeds the predefined region radius, given by HdVR, the point is accepted; otherwise, a new trial point is generated. When one of the best obtained trial solutions is sufficiently close to a global minimum, or its value has not been changed for a specified number of iterations, the DTS algorithm leaves the exploration and diversification search procedures and invokes an intensification procedure. Here, the HJ local search is used to compute a solution still closer to the global minimum (see the next section).


## 3.6. Hooke-and-Jeeves Local Search

A pattern search method directs the search towards a minimizer using a pattern of specific number points. At leastn+1 points must be provided by the pattern, where n is the number of variables. Based on a current approximation, the HJ method uses a pattern search strategy to define a trial solution that gives an improvement in the objective function value, when compared with the current point [25]. Let the current approximation be denoted by p(ki), where ki represents the iteration counter in this iterative process. Then, a trial point, vi, is generated along a search direction (starting from the current point) with a certain step size Δ(ki)>0 as follows:(11)vi=pki+Δkidi,where di is the search direction chosen from a finite set D of positive spanning directions in Rn. The most used set contains the 2n coordinate directions, defined as the positive and negative unit coordinate vectors D=e1,…,en,-e1,…,-en. The most important property of this set is that at least one of the coordinate directions is a descent direction for the objective J, when p(ki) is not a stationary point of J. When the search fails to generate a trial point that is better than p(ki), the iteration is calledunsuccessful, the step size Δ(ki) is halved so that a refined search can be carried out, and p(ki+1)←p(ki). If Δ(ki) falls below a given stopping tolerance, the algorithm terminates and the current p(ki) is considered an approximate minimizer of J. However, when at the end of each iteration the objective function value has reduced, then the iteration is calledsuccessful, Δ(ki) is not changed, and p(ki+1)←vi.For a fair comparison, and to improve the quality of the produced solutions, we propose the implementation of the HJ intensification phase with all the above-mentioned metaheuristics. Based on the final solution provided by the metaheuristic, the HJ local search algorithm is invoked and allowed to run for5n iterations. The number of function evaluations required by the HJ local search is added to the function evaluations of the metaheuristic, herein denoted by NFEval, to produce the total number of evaluations of the run.


## 4. Numerical Experiments

This section aims to present and analyze the statistical significance of the numerical results that were obtained when the metaheuristics are used to solve the DMbPE problems in the sequential direct method context. The dynamic system models were coded in the Matlab programming language and the computational application of the sequential direct method with the FA, HS, DE, ABC, and DTS codes was developed in Matlab programming environment. The computational tests were performed on a PC with a 2.2 GHz Core i7-2670QM and 8 GB of RAM. The parameter values for the tested algorithms are set as shown in Table6. Other values have been tested, in particular for the parameter NP in the FA, DE, and ABC and the parameter HMS in HS, but those used seem adequate with reduced computational effort.Table 6
Parameter values.
FA
N
P
min
⁡
{
5
n
,
50
}
α (initial)
0.25
γ (constant)
1
β
0
1
β
m
i
n
0.2
HS
H
M
S
min
⁡
{
2
n
,
10
}
H
M
C
R
0.95
DE
N
P
min
⁡
4
n
,
50
F, CR
0.9
Mutation
DE/rand-to-best/1
ABC
N
P
min
⁡
{
4
n
,
50
}
“limit”
100
DTS
H
1
f
e
d
g
e
0.1
r
d
V
R
2Nine case study problems have been selected to analyze the performance of the described stochastic algorithms. A comparison is made considering the quality of the solutions and the time spent to reach the solution after a threshold number of function evaluations. The problems and the experimental data can be found in the Appendix.Table7 shows Javg (the average of the obtained final solutions), St.D. (the standard deviations of the solutions), and Tavg (the average time in seconds) after running each instance (inst.) 10 times with each metaheuristic (MH). For these experiments, each algorithm is terminated solely with the condition NFEval.>NFmax, where we have tested three values of NFmax: 50n, 100n, and 200n. The goal here is to analyze the quality of the obtained solutions. Since the target objective value is 0, the lower the value of Javg the better.Table 7
Comparison ofJavg, St.D., and Tavg.
Inst.
n
MH
N
F
m
a
x
=
50
n
N
F
m
a
x
=
100
n
N
F
m
a
x
=
200
n
J
a
v
g
St.D.
T
a
v
g
J
a
v
g
St.D.
T
a
v
g
J
a
v
g
St.D.
T
a
v
g
α
-PIN
5
FA
1.99583
E
+
01
9.98
E
-
02
2.09
E
+
01
1.99145
E
+
01
9.91
E
-
02
3.27
E
+
01
2.00343
E
+
01
4.28
E
-
01
5.45
E
+
01
HS
1.99172
E
+
01
1.20
E
-
01
2.31
E
+
01
1.98832
E
+
01
1.50
E
-
02
3.16
E
+
01
1.98974
E
+
01
5.27
E
-
02
6.01
E
+
01
DE
2.02428
E
+
01
4.98
E
-
01
2.51
E
+
01
2.03504
E
+
01
1.10
E
+
00
3.20
E
+
01
1.98827
E
+
01
1.91
E
-
02
5.30
E
+
01
ABC
4.40554
E
+
02
6.38
E
+
02
2.02
E
+
01
8.97875
E
+
01
7.64
E
+
01
3.08
E
+
01
2.30471
E
+
01
3.15
E
+
00
5.76
E
+
01
DTS
2.02791
E
+
01
8.23
E
-
01
2.56
E
+
01
2.00559
E
+
01
3.93
E
-
01
3.94
E
+
01
2.00243
E
+
01
3.56
E
-
01
1.03
E
+
02
BEL_A
2
FA
2.21943
E
+
01
2.02
E
-
02
4.01
E
+
00
2.21867
E
+
01
6.26
E
-
03
6.34
E
+
00
2.21922
E
+
01
2.13
E
-
02
1.08
E
+
01
HS
2.22921
E
+
01
2.73
E
-
01
4.18
E
+
00
2.22393
E
+
01
1.68
E
-
01
6.45
E
+
00
2.21886
E
+
01
9.73
E
-
03
1.13
E
+
01
DE
2.22109
E
+
01
3.72
E
-
02
4.00
E
+
00
2.21871
E
+
01
6.95
E
-
03
6.18
E
+
00
2.21815
E
+
01
1.56
E
-
04
1.08
E
+
01
ABC
4.56940
E
+
01
2.83
E
+
01
4.06
E
+
00
2.40930
E
+
01
2.86
E
+
00
6.50
E
+
00
2.23062
E
+
01
2.39
E
-
01
1.12
E
+
01
DTS
2.22325
E
+
01
7.77
E
-
02
5.51
E
+
00
2.23786
E
+
01
4.51
E
-
01
9.47
E
+
00
2.22101
E
+
01
2.84
E
-
02
1.81
E
+
01
BEL_B
2
FA
2.21928
E
+
01
1.31
E
-
02
3.67
E
+
00
2.21867
E
+
01
7.67
E
-
03
6.43
E
+
00
2.21821
E
+
01
3.58
E
-
04
1.09
E
+
01
HS
2.22161
E
+
01
4.15
E
-
02
3.67
E
+
00
2.22001
E
+
01
2.55
E
-
02
6.02
E
+
00
2.21870
E
+
01
6.57
E
-
03
1.12
E
+
01
DE
2.21918
E
+
01
1.54
E
-
02
3.80
E
+
00
2.21824
E
+
01
1.52
E
-
03
6.12
E
+
00
2.21814
E
+
01
3.01
E
-
07
1.09
E
+
01
ABC
3.13748
E
+
01
1.32
E
+
01
3.54
E
+
00
2.49722
E
+
01
6.36
E
+
00
6.10
E
+
00
2.25338
E
+
01
1.11
E
+
00
1.11
E
+
01
DTS
2.22010
E
+
01
1.44
E
-
02
4.16
E
+
00
2.22058
E
+
01
2.59
E
-
02
7.01
E
+
00
2.22025
E
+
01
3.51
E
-
02
1.23
E
+
01
CAT
3
FA
2.6559
E
-
03
3.59
E
-
07
1.07
E
+
01
2.6558
E
-
03
1.52
E
-
07
1.73
E
+
01
2.6557
E
-
03
1.25
E
-
08
4.85
E
+
01
HS
2.6607
E
-
03
1.50
E
-
05
1.03
E
+
01
2.6558
E
-
03
1.31
E
-
07
1.72
E
+
01
2.6558
E
-
03
2.17
E
-
07
3.38
E
+
01
DE
2.6566
E
-
03
1.98
E
-
06
1.10
E
+
01
2.6557
E
-
03
1.46
E
-
07
1.68
E
+
01
2.6557
E
-
03
6.03
E
-
09
3.29
E
+
01
ABC
1.1355
E
-
02
8.65
E
-
03
9.73
E
+
00
4.2229
E
-
03
1.65
E
-
03
1.65
E
+
01
3.0625
E
-
03
4.94
E
-
04
3.46
E
+
01
DTS
2.6578
E
-
03
4.72
E
-
06
1.12
E
+
01
2.6560
E
-
03
5.19
E
-
07
1.78
E
+
01
2.6570
E
-
03
2.55
E
-
06
3.46
E
+
01
IRR1
2
FA
2.9302
E
-
06
1.83
E
-
06
4.30
E
+
00
2.1233
E
-
06
1.38
E
-
06
6.42
E
+
00
1.3424
E
-
06
3.49
E
-
07
1.20
E
+
01
HS
3.9736
E
-
06
3.22
E
-
06
4.11
E
+
00
2.2205
E
-
06
1.21
E
-
06
6.49
E
+
00
1.4775
E
-
06
3.14
E
-
07
1.17
E
+
01
DE
4.3658
E
-
06
6.02
E
-
06
4.21
E
+
00
1.5015
E
-
06
5.28
E
-
07
6.54
E
+
00
1.1871
E
-
06
2.86
E
-
09
1.18
E
+
01
ABC
8.4434
E
-
02
1.70
E
-
01
3.81
E
+
00
2.3661
E
-
04
4.22
E
-
04
6.76
E
+
00
2.3818
E
-
03
7.53
E
-
03
1.21
E
+
01
DTS
7.2443
E
-
06
8.14
E
-
06
4.06
E
+
00
2.5000
E
-
06
1.24
E
-
06
6.67
E
+
00
5.4046
E
-
06
5.09
E
-
06
1.20
E
+
01
IRR2
2
FA
6.9479
E
-
05
7.39
E
-
05
5.27
E
+
00
2.5916
E
-
05
3.78
E
-
05
8.63
E
+
00
6.5898
E
-
06
2.81
E
-
06
1.56
E
+
01
HS
2.6359
E
-
04
5.05
E
-
04
5.20
E
+
00
3.9423
E
-
05
5.12
E
-
05
8.50
E
+
00
2.6067
E
-
05
2.47
E
-
05
1.50
E
+
01
DE
2.0207
E
-
05
4.06
E
-
05
5.37
E
+
00
6.6501
E
-
06
2.66
E
-
06
8.64
E
+
00
4.4153
E
-
06
7.19
E
-
09
1.55
E
+
01
ABC
6.5399
E
-
03
1.50
E
-
02
4.94
E
+
00
1.1674
E
-
03
3.61
E
-
03
8.59
E
+
00
4.5931
E
-
06
5.31
E
-
07
1.55
E
+
01
DTS
4.5181
E
-
05
5.14
E
-
05
4.87
E
+
00
4.5081
E
-
05
3.55
E
-
05
7.97
E
+
00
5.6111
E
-
05
8.44
E
-
05
1.42
E
+
01
REV_A
4
FA
8.9843
E
-
06
2.48
E
-
05
1.28
E
+
01
1.2775
E
-
06
1.44
E
-
06
1.90
E
+
01
2.1194
E
-
06
3.47
E
-
06
3.17
E
+
01
HS
3.6938
E
-
05
9.88
E
-
05
1.29
E
+
01
4.3732
E
-
05
1.08
E
-
04
1.92
E
+
01
1.4107
E
-
06
1.41
E
-
06
3.17
E
+
01
DE
5.5673
E
-
07
3.03
E
-
07
3.09
E
+
01
1.5841
E
-
06
2.92
E
-
06
1.94
E
+
01
6.1260
E
-
07
5.65
E
-
07
3.00
E
+
01
ABC
4.6219
E
-
02
3.96
E
-
02
9.18
E
+
00
1.9570
E
-
02
4.61
E
-
02
1.68
E
+
01
4.0868
E
-
03
2.59
E
-
03
3.01
E
+
01
DTS
5.8993
E
-
04
1.10
E
-
03
1.28
E
+
01
7.9723
E
-
05
7.71
E
-
05
1.95
E
+
01
3.5410
E
-
05
4.42
E
-
05
3.24
E
+
01
REV_B
4
FA
1.6175
E
-
03
8.04
E
-
05
1.27
E
+
01
1.5897
E
-
03
5.59
E
-
06
1.86
E
+
01
1.5884
E
-
03
1.39
E
-
06
3.15
E
+
01
HS
1.6440
E
-
03
1.07
E
-
04
1.28
E
+
01
1.5980
E
-
03
1.01
E
-
05
1.94
E
+
01
1.5884
E
-
03
8.81
E
-
07
3.17
E
+
01
DE
1.5877
E
-
03
3.90
E
-
07
3.07
E
+
01
1.5879
E
-
03
7.44
E
-
07
2.95
E
+
01
1.5878
E
-
03
5.17
E
-
07
2.98
E
+
01
ABC
3.7122
E
-
02
4.25
E
-
02
9.11
E
+
00
1.5726
E
-
02
1.16
E
-
02
1.67
E
+
01
7.1333
E
-
03
6.73
E
-
03
3.06
E
+
01
DTS
1.8431
E
-
03
3.23
E
-
04
1.29
E
+
01
1.6313
E
-
03
6.82
E
-
05
1.93
E
+
01
1.7760
E
-
03
2.33
E
-
04
3.27
E
+
01
THEO1
3
FA
1.0938
E
-
03
8.53
E
-
06
9.47
E
+
00
1.0842
E
-
03
8.85
E
-
06
1.47
E
+
01
1.0958
E
-
03
9.80
E
-
06
2.51
E
+
01
HS
1.0900
E
-
03
1.61
E
-
05
9.11
E
+
00
1.0842
E
-
03
1.61
E
-
05
1.42
E
+
01
1.0861
E
-
03
1.55
E
-
05
2.46
E
+
01
DE
1.0706
E
-
03
1.42
E
-
08
8.64
E
+
00
1.0764
E
-
03
9.19
E
-
06
1.35
E
+
01
1.0763
E
-
03
9.19
E
-
06
2.33
E
+
01
ABC
1.0934
E
-
03
1.49
E
-
05
8.07
E
+
00
1.0809
E
-
03
9.90
E
-
06
1.38
E
+
01
1.0766
E
-
03
7.55
E
-
06
2.41
E
+
01
DTS
1.0978
E
-
03
1.05
E
-
05
1.05
E
+
01
1.0860
E
-
03
8.10
E
-
06
1.67
E
+
01
1.0841
E
-
03
1.31
E
-
05
2.88
E
+
01
THEO2
2
FA
3.5953
E
-
02
1.52
E
-
04
6.38
E
+
00
3.5800
E
-
02
1.38
E
-
05
1.08
E
+
01
3.5784
E
-
02
4.51
E
-
06
1.88
E
+
01
HS
3.5976
E
-
02
1.92
E
-
04
6.98
E
+
00
3.5866
E
-
02
9.75
E
-
05
1.14
E
+
01
3.5809
E
-
02
4.50
E
-
05
2.03
E
+
01
DE
3.5829
E
-
02
6.46
E
-
05
6.74
E
+
00
3.5782
E
-
02
8.44
E
-
07
1.04
E
+
01
3.5809
E
-
02
4.50
E
-
05
2.03
E
+
01
ABC
4.6244
E
-
02
2.17
E
-
02
8.47
E
+
00
2.1159
E
-
01
5.52
E
-
01
1.30
E
+
01
4.2209
E
-
02
2.03
E
-
02
2.12
E
+
01
DTS
3.5863
E
-
02
1.24
E
-
04
7.92
E
+
00
3.5876
E
-
02
1.15
E
-
04
1.39
E
+
01
3.5833
E
-
02
5.42
E
-
05
2.56
E
+
01
VOL_A
2
FA
2.1700
E
-
01
2.52
E
-
01
3.13
E
+
01
1.6923
E
-
01
2.28
E
-
01
4.50
E
+
01
1.5620
E
-
01
2.28
E
-
01
8.31
E
+
01
HS
1.0265
E
-
01
1.66
E
-
01
2.90
E
+
01
9.8232
E
-
02
1.61
E
-
01
5.51
E
+
01
3.4138
E
-
02
2.19
E
-
02
8.89
E
+
01
DE
2.0002
E
-
02
1.08
E
-
02
3.31
E
+
01
1.6915
E
-
02
9.57
E
-
03
4.76
E
+
01
6.2283
E
-
02
1.48
E
-
01
9.26
E
+
01
ABC
1.7601
E
-
01
2.94
E
-
01
2.73
E
+
01
1.0289
E
-
01
1.56
E
-
01
4.53
E
+
01
3.8400
E
-
02
3.45
E
-
02
8.02
E
+
01
DTS
2.9206
E
-
02
1.73
E
-
02
2.96
E
+
01
8.1213
E
-
02
1.46
E
-
01
4.70
E
+
01
2.5780
E
-
02
6.98
E
-
03
8.79
E
+
01
VOL_B
2
FA
1.0574
E
+
00
1.77
E
-
01
3.35
E
+
01
1.0277
E
+
00
1.66
E
-
01
7.16
E
+
01
9.7332
E
-
01
1.61
E
-
01
6.59
E
+
01
HS
9.5826
E
-
01
1.03
E
-
01
2.75
E
+
01
1.0034
E
+
00
1.58
E
-
01
4.18
E
+
01
9.7446
E
-
01
1.33
E
-
01
7.41
E
+
01
DE
9.3036
E
-
01
3.68
E
-
02
3.26
E
+
01
9.0152
E
-
01
2.33
E
-
02
3.68
E
+
01
9.3849
E
-
01
1.11
E
-
01
7.89
E
+
01
ABC
1.0558
E
+
00
2.30
E
-
01
2.54
E
+
01
9.7174
E
-
01
1.55
E
-
01
4.56
E
+
01
9.1724
E
-
01
5.10
E
-
02
6.42
E
+
01
DTS
9.7639
E
-
01
1.14
E
-
01
2.81
E
+
01
9.7674
E
-
01
1.04
E
-
01
4.54
E
+
01
8.9907
E
-
01
2.99
E
-
02
7.78
E
+
01To analyze the statistical significance of the results, we use the Matlab functionfriedman that performs a nonparametric statistical test for multiple comparisons. The Friedman test aims to determine significant differences in the mean for one independent variable, for instance, Javg or Tavg, with different levels that correspond to the five metaheuristics, and a dependent variable (corresponding to matched groups, herein taken as the 12 instances) [35]. The null hypothesis in this test is that the mean ranks assigned to the results of the metaheuristics under testing are not statistically different. Let k be the number of metaheuristics to be ranked and N be the number of groups. We note here that the distribution of the Friedmanstatistic, here denoted by QF, approaches the ordinary χ2 distribution as N increases. The exact distribution of QF has been obtained for special cases of k and N (the reader is referred to [35] and the references cited in [36]). It has been concluded in [35] that “when the number of groups is moderately large (say greater than 5 for four or more ranks) the significance of the Friedmanstatistic can be tested by reference to the available χ2 tables.” However, in [36], it is argued that “the usual χ2 approximation is grossly inaccurate for most commonly used combinations of k and N.” Thus, hereafter, we use anotherstatistic that has been recommended in [36, 37] and depends on QF, as well as on both k and N:(12)FF=N-1QFNk-1-QF,which is distributed according to the F-distribution with k-1 and (k-1)(N-1) degrees of freedom.When applied toJavg for NFmax=50n, QF=25.27 and the value of FF is 12.229. With k=5 and N=12, FF is distributed according to the F-distribution with 4 and 44 degrees of freedom. For a significance level of 5%, the critical value F(4,44)0.05 is 2.59 (computed by linear interpolation between F(4,40)0.05 and F(4,50)0.05). Since 12.229>F(4,44)0.05, we have enough evidence to reject the null hypothesis of “no significant differences on mean ranks”; that is, the observed differences between the five distributions of Javg values are statistically significant. Pairwise comparisons may be carried out to determine which mean ranks are significantly different. The Matlab function multcompare is applied. The estimates of the 95% confidence intervals are shown in Figure 1(a) for each case under testing. The distributions of Javg are significantly different if their intervals are disjoint and are not significantly different if their intervals overlap. Hence, we conclude that the mean rank produced by DE is significantly different from that of ABC. For the other pairs of comparisons, there are no significant differences on the mean ranks.Figure 1
Estimates of the 95% confidence intervals forJavg and Tavg, using NFmax⁡=50n.
(a)
J
a
v
g metric
(b)
T
a
v
g metricWhen the test is applied toTavg for NFmax=50n, QF=14.97 and FF=4.985 which exceeds the value of F(4,44)0.05, indicating that the null hypothesis is rejected at a significance level of 5%, and we conclude that the distributions of Tavg values have statistically significant differences. Figure 1(b) shows the estimates of the 95% confidence intervals. We conclude that the mean rank of ABC is significantly different from those of DE and DTS.A similar analysis is made forJavg and Tavg, when NFmax=100n. We obtain for the statistics QF=25.82 with FF=12.805 and QF=7.47 with FF=2.027, respectively. Hence, we may conclude that there is enough evidence to reject the null hypothesis of “no significant differences on mean ranks” of the distributions of Javg, at a significance level of 5%. From Figure 2(a), we conclude that the mean ranks of (the distribution of Javg values) ABC and DTS are significantly different from the mean ranks of DE. For the Tavg distributions, we conclude that there is not enough evidence to reject the hypothesis of “no significant differences on mean ranks,” at a significance level of 5%. Figure 2(b) shows the corresponding estimates of the 95% confidence intervals and it is noticed that all pairs have their intervals overlapping.Figure 2
Estimates of the 95% confidence intervals forJavg and Tavg, using NFmax⁡=100n.
(a)
J
a
v
g metric
(b)
T
a
v
g metricWhen the statistical analysis is extended to the distributions ofJavg and Tavg values, for NFmax=200n, the statistic QF gives 16.34 (with FF=5.677) and 10.64 (with FF=3.133), respectively. Hence, we have evidence to reject the null hypothesis relative to the Javg values, at a significance level of 5%. Figure 3(a) shows the corresponding estimates of the 95% confidence intervals (ABC and DTS have mean ranks significantly different from DE). As far as the distributions of Tavg values are concerned, we conclude that there is enough evidence to reject the hypothesis of “no significant differences on mean ranks” at a significance level of 5%. From Figure 3(b), we conclude that the differences are on the mean ranks of DE and DTS. We note that, at a significance level of 1%, a smaller probability of making a wrong decision, where the critical value F(4,44)0.01 is 3.79 (computed by linear interpolation between F(4,40)0.01 and F(4,50)0.01), there is no evidence to reject the null hypothesis.Figure 3
Estimates of the 95% confidence intervals forJavg and Tavg, using NFmax=200n.
(a)
J
a
v
g metric
(b)
T
a
v
g metricWe now show in Table8 the best solution obtained by each metaheuristic, Jbest (after the 10 runs), and the time required to reach that value, Tbest, for the three values of NFmax. Based on the metrics Jbest and Tbest, we use a graphical procedure to visualize the performance differences among the results produced by the five metaheuristics, in relative terms on the 12 instances, known as performance profiles [38]. Each plot reports (on the vertical axis) the percentage of problems solved with each metaheuristic that is within a certain threshold, τ (on the horizontal axis), of the best result. The performance profiles of the five distributions of Jbest values for NFmax=50n are shown in Figure 4(a). Figure 4(b) corresponds to the metric Tbest.Table 8
Best results (Tbest and Jbest).
Inst.
MH
N
F
m
a
x
=
50
n
N
F
m
a
x
=
100
n
N
F
m
a
x
=
200
n
T
b
e
s
t
J
b
e
s
t
T
b
e
s
t
J
b
e
s
t
T
b
e
s
t
J
b
e
s
t
α
-PIN
FA
2.17
E
+
01
1.987279
E
+
01
3.08
E
+
01
1.987226
E
+
01
4.80
E
+
01
1.987264
E
+
01
HS
2.15
E
+
01
1.987260
E
+
01
3.04
E
+
01
1.987228
E
+
01
5.18
E
+
01
1.987230
E
+
01
DE
2.29
E
+
01
1.987271
E
+
01
2.92
E
+
01
1.987296
E
+
01
5.18
E
+
01
1.987263
E
+
01
ABC
1.51
E
+
01
7.990355
E
+
01
2.54
E
+
01
2.567832
E
+
01
4.82
E
+
01
2.036694
E
+
01
DTS
2.31
E
+
01
1.987234
E
+
01
3.83
E
+
01
1.987266
E
+
01
6.91
E
+
01
1.987274
E
+
01
BEL_A
FA
4.29
E
+
00
2.218250
E
+
01
5.97
E
+
00
2.218155
E
+
01
1.10
E
+
01
2.218149
E
+
01
HS
3.95
E
+
00
2.218227
E
+
01
6.15
E
+
00
2.218149
E
+
01
1.14
E
+
01
2.218311
E
+
01
DE
4.48
E
+
00
2.218143
E
+
01
6.05
E
+
00
2.218147
E
+
01
1.06
E
+
01
2.218142
E
+
01
ABC
3.44
E
+
00
2.227154
E
+
01
5.78
E
+
00
2.218221
E
+
01
1.07
E
+
01
2.218150
E
+
01
DTS
4.59
E
+
00
2.218762
E
+
01
9.52
E
+
00
2.218738
E
+
01
1.76
E
+
01
2.218458
E
+
01
BEL_B
FA
5.09
E
+
00
2.218269
E
+
01
6.07
E
+
00
2.218152
E
+
01
1.09
E
+
01
2.218170
E
+
01
HS
3.52
E
+
00
2.218164
E
+
01
6.04
E
+
00
2.218142
E
+
01
1.12
E
+
01
2.218159
E
+
01
DE
3.68
E
+
00
2.218181
E
+
01
6.26
E
+
00
2.218142
E
+
01
1.11
E
+
01
2.218142
E
+
01
ABC
3.17
E
+
00
2.218809
E
+
01
5.89
E
+
00
2.218259
E
+
01
1.15
E
+
01
2.218151
E
+
01
DTS
4.32
E
+
00
2.218208
E
+
01
7.33
E
+
00
2.218152
E
+
01
1.23
E
+
01
2.218187
E
+
01
CAT
FA
1.03
E
+
01
2.655668
E
-
03
1.76
E
+
01
2.655668
E
-
03
4.78
E
+
01
2.655668
E
-
03
HS
1.01
E
+
01
2.655686
E
-
03
1.61
E
+
01
2.655669
E
-
03
3.23
E
+
01
2.655669
E
-
03
DE
1.06
E
+
01
2.655683
E
-
03
1.62
E
+
01
2.655667
E
-
03
3.11
E
+
01
2.655667
E
-
03
ABC
9.82
E
+
00
2.666711
E
-
03
1.65
E
+
01
2.946387
E
-
03
3.45
E
+
01
2.662231
E
-
03
DTS
1.05
E
+
01
2.655667
E
-
03
1.73
E
+
01
2.655680
E
-
03
3.40
E
+
01
2.655672
E
-
03
IRR1
FA
4.45
E
+
00
1.221743
E
-
06
6.51
E
+
00
1.228962
E
-
06
1.18
E
+
01
1.187445
E
-
06
HS
4.29
E
+
00
1.375865
E
-
06
6.53
E
+
00
1.190046
E
-
06
1.15
E
+
01
1.193409
E
-
06
DE
4.49
E
+
00
1.188136
E
-
06
6.41
E
+
00
1.194300
E
-
06
1.18
E
+
01
1.185902
E
-
06
ABC
4.33
E
+
00
4.204506
E
-
05
6.90
E
+
00
1.271379
E
-
06
1.21
E
+
01
1.185851
E
-
06
DTS
4.22
E
+
00
1.333455
E
-
06
6.77
E
+
00
1.228666
E
-
06
1.18
E
+
01
1.228072
E
-
06
IRR2
FA
5.49
E
+
00
5.793696
E
-
06
8.67
E
+
00
5.290805
E
-
06
1.55
E
+
01
4.532302
E
-
06
HS
5.40
E
+
00
7.799776
E
-
06
8.10
E
+
00
7.581183
E
-
06
1.54
E
+
01
4.946930
E
-
06
DE
5.65
E
+
00
4.951955
E
-
06
8.65
E
+
00
4.512205
E
-
06
1.53
E
+
01
4.411032
E
-
06
ABC
5.15
E
+
00
1.539352
E
-
05
8.64
E
+
00
4.473528
E
-
06
1.52
E
+
01
4.411030
E
-
06
DTS
4.82
E
+
00
5.176208
E
-
06
8.00
E
+
00
8.305471
E
-
06
1.41
E
+
01
5.550882
E
-
06
REV_A
FA
1.30
E
+
01
2.390974
E
-
07
1.95
E
+
01
2.301914
E
-
07
3.13
E
+
01
1.933663
E
-
07
HS
1.29
E
+
01
4.456872
E
-
07
1.95
E
+
01
2.772230
E
-
07
3.09
E
+
01
2.864920
E
-
07
DE
3.13
E
+
01
1.960909
E
-
07
1.78
E
+
01
2.060756
E
-
07
3.03
E
+
01
2.034098
E
-
07
ABC
1.03
E
+
01
1.184856
E
-
02
1.64
E
+
01
3.235180
E
-
04
3.05
E
+
01
4.414293
E
-
04
DTS
1.28
E
+
01
3.003054
E
-
07
1.90
E
+
01
2.925082
E
-
07
3.22
E
+
01
4.346120
E
-
07
REV_B
FA
1.16
E
+
01
1.587466
E
-
03
1.70
E
+
01
1.587489
E
-
03
3.22
E
+
01
1.587486
E
-
03
HS
1.30
E
+
01
1.587876
E
-
03
1.86
E
+
01
1.587498
E
-
03
3.08
E
+
01
1.587476
E
-
03
DE
3.07
E
+
01
1.587467
E
-
03
2.88
E
+
01
1.587501
E
-
03
2.89
E
+
01
1.587516
E
-
03
ABC
7.27
E
+
00
7.753336
E
-
03
1.70
E
+
01
2.810094
E
-
03
3.03
E
+
01
1.816843
E
-
03
DTS
1.32
E
+
01
1.587997
E
-
03
1.90
E
+
01
1.587864
E
-
03
3.36
E
+
01
1.587618
E
-
03
THEO1
FA
9.28
E
+
00
1.089660
E
-
03
1.35
E
+
01
1.070636
E
-
03
2.58
E
+
01
1.089664
E
-
03
HS
8.95
E
+
00
1.070635
E
-
03
1.39
E
+
01
1.070637
E
-
03
2.33
E
+
01
1.070635
E
-
03
DE
8.21
E
+
00
1.070638
E
-
03
1.31
E
+
01
1.070635
E
-
03
2.25
E
+
01
1.070635
E
-
03
ABC
7.15
E
+
00
1.071921
E
-
03
1.23
E
+
01
1.071302
E
-
03
2.30
E
+
01
1.070637
E
-
03
DTS
1.03
E
+
01
1.089659
E
-
03
1.61
E
+
01
1.070635
E
-
03
2.78
E
+
01
1.070638
E
-
03
THEO2
FA
6.22
E
+
00
3.580257
E
-
02
1.10
E
+
01
3.578306
E
-
02
1.89
E
+
01
3.578128
E
-
02
HS
6.78
E
+
00
3.579577
E
-
02
1.13
E
+
01
3.578748
E
-
02
2.03
E
+
01
3.578180
E
-
02
DE
6.45
E
+
00
3.578201
E
-
02
1.02
E
+
01
3.578109
E
-
02
2.03
E
+
01
3.578180
E
-
02
ABC
7.99
E
+
00
3.608823
E
-
02
1.24
E
+
01
3.580068
E
-
02
2.09
E
+
01
3.578108
E
-
02
DTS
8.76
E
+
00
3.578341
E
-
02
1.42
E
+
01
3.578354
E
-
02
2.54
E
+
01
3.578140
E
-
02
VOL_A
FA
1.37
E
+
01
7.404584
E
-
03
2.56
E
+
01
9.624342
E
-
03
3.33
E
+
01
3.602099
E
-
03
HS
1.60
E
+
01
1.550778
E
-
02
3.04
E
+
01
2.731927
E
-
02
3.82
E
+
01
1.388033
E
-
02
DE
1.19
E
+
01
3.787360
E
-
03
1.81
E
+
01
3.574318
E
-
03
4.06
E
+
01
3.574082
E
-
03
ABC
2.37
E
+
01
5.221512
E
-
03
2.73
E
+
01
4.371459
E
-
03
4.51
E
+
01
4.326559
E
-
03
DTS
2.33
E
+
01
1.018590
E
-
02
4.11
E
+
01
5.239623
E
-
03
8.39
E
+
01
1.922888
E
-
02
VOL_B
FA
1.36
E
+
01
9.124563
E
-
01
3.72
E
+
01
8.667088
E
-
01
4.07
E
+
01
8.684604
E
-
01
HS
1.66
E
+
01
8.700917
E
-
01
2.48
E
+
01
8.810708
E
-
01
3.90
E
+
01
8.697327
E
-
01
DE
1.54
E
+
01
8.709003
E
-
01
1.89
E
+
01
8.703102
E
-
01
4.67
E
+
01
8.667010
E
-
01
ABC
1.51
E
+
01
8.955672
E
-
01
2.02
E
+
01
8.671336
E
-
01
6.42
E
+
01
8.700530
E
-
01
DTS
2.66
E
+
01
8.919989
E
-
01
3.87
E
+
01
9.019120
E
-
01
7.12
E
+
01
8.679623
E
-
01Figure 4
Performance profiles for the five metaheuristics, usingNFmax=50n.
(a)
J
best metric
(b)
T
best metricThe higher the percentage, the better. A higher value forτ=1 means that the metaheuristic achieves the lowest Jbest value (or the smallest Tbest) mostly.Therefore, the metaheuristic DE is the most successful since it has the highest probability of achieving the lowestJ value, when NFmax=50n. From the plot, the probability that DE wins on a given instance is 0.5 since it produces the lowest Jbest in 6 of the 12 instances. It is followed by HS that has a probability of winning (prob_w) of 0.25, by DTS with prob_w = 0.17, and by FA with prob_w = 0.08. We note that the differences on the produced Jbest are so small that the plots of the cumulative distribution are packed around τ=1. From Figure 4(b), it is possible to conclude that the fastest metaheuristic is ABC since the probability of achieving the smallest value of T (time in seconds) is 0.58.Similarly, Figures5(a) and 5(b) contain the performance profiles of the five distributions of Jbest and Tbest, respectively, when NFmax=100n. From the figure on the left, we may conclude that DE is the most successful since it has prob_w = 0.58, followed by FA with prob_w = 0.25 (achieving the best results in 3 out of 12 instances), by HS with prob_w = 0.17, and by ABC and DTS with prob_w = 0.08. As far as Tbest is concerned, it is possible to conclude that the metaheuristic with the highest probability of requiring the lowest time is ABC followed by DE.Figure 5
Performance profiles for the five metaheuristics, usingNFmax=100n.
(a)
J
best metric
(b)
T
best metricFinally, Figures6(a) and 6(b) display the performance profiles of the metrics Jbest and Tbest, respectively, when NFmax=200n. From these comparisons, we may conclude that DE is the most successful in reaching the lowest Jbest, with prob_w = 0.5, followed by HS and ABC with prob_w = 0.25 and finally the metaheuristics FA and DTS with prob_w = 0.08. When the profiles of Tbest are analyzed, it is possible to conclude that the metaheuristic DE is in general faster than the others but for some problems is outstripped by FA and HS.Figure 6
Performance profiles for the five metaheuristics, usingNFmax=200n.
(a)
J
best metric
(b)
T
best metricTables9–19 contain comparative results relative to 11 instances. (The results obtained for instance VOL_A are not used in the comparison since the experimental data have been herein generated for this study.) Tables 9–19 show the best and the average J values (computed from the 10 runs), as well as the average CPU time and function evaluations, Tavg and NFEavg, respectively. Our results always appear along the first row of each table and are taken from the metaheuristic that produced the lowest Jbest with the least function evaluations. This is identified in the table caption. (We note that the implemented metaheuristics have been enhanced with the HJ local search, and the reported NFEavg considers the sum of NFmax with the function evaluations required by the HJ search.) In the tables, “–” means that the information is not available in the cited paper and “n.a.” means “not applicable.” The comparisons use the results reported in the papers [2, 4, 10, 12–14, 16–19, 26]. Results in [16] are based on a trigonometric version of the DE, and in [17], a modified version of the DE is implemented. The authors in [12] use a modified version of the metaheuristic scatter search (SSm) with a local search as an improvement method. The results shown in the table use the function fmincon from Matlab. SSm uses a population of 100 points and a penalty technique with a weight of 1000 to penalize infeasible solutions. In [19], the sequential approach is implemented with an interval-based method for the NLP optimization. The therein obtained solution is an ϵ-global solution. Global solutions are produced by the method of orthogonal collocation on finite elements presented in [26]. This strategy is implemented within Matlab code Dynopt after 100 runs of a multistart approach. (We note that a direct comparison of the CPU time cannot be done since the computational platforms are different.)Table 9
Comparative results forα-PIN (FA/100n).
Source
J
b
e
s
t
J
a
v
g
T
a
v
g
N
F
E
a
v
g
1.987226
E
+
01
1.991449
E
+
01
3.27
E
+
01
753
[2]♭
1.987
E
+
01
–
1.22
E
+
02
(9518 iter.)
[12]†
1.9872
E
+
01
2.4747
E
+
01
0.41
E
+
02
1163
[13]§
1.98772
E
+
01
2.56777
E
+
01
–
9148
[14]‡
1.987
E
+
01
n.a.
8.92
E
+
03
(2 iter.)
♭PC: Pentium 4, 1.80 GHz. †PC: Pentium 4, 3.06 GHz. §PC with a 2.7 GHz Core i7-4600U and 8 GB of memory. ‡PC: AMD Athlon II, 2.99 GHz.Table 10
Comparative results forBEL_A (loose bounds on parameters) (DE/200n).
Source
J
b
e
s
t
J
a
v
g
T
a
v
g
N
F
E
a
v
g
2.218142
E
+
01
2.218148
E
+
01
1.08
E
+
01
434
[10]†
2.25684
E
+
01
n.a.
9.71
E
+
03
(10000 iter.)
†Result with relative convergence tolerance of 7.5E-01. Computer HP J2240.Table 11
Comparative results forBEL_B (tight bounds on parameters) (HS/100n).
Source
J
b
e
s
t
J
a
v
g
T
a
v
g
N
F
E
a
v
g
2.218142
E
+
01
2.220014
E
+
01
6.02
E
+
00
245
[10]†
2.25684
E
+
01
n.a.
7.84
E
+
03
(9028 iter.)
†Result with relative convergence tolerance of 1.0E-03. Computer HP J2240.Table 12
Comparative results forCAT (DTS/50n).
Source
J
b
e
s
t
J
a
v
g
T
a
v
g
N
F
E
a
v
g
2.655667
E
-
03
2.657778
E
-
03
1.12
E
+
01
256
[4]
2.65567
E
-
03
n.a.
–
–
[10]♭
2.6384
E
-
03
n.a.
1.58
E
+
02
(26 iter.)
[16]§
2.655666
E
-
03
2.655666
E
-
03
1.81
E
+
00
–
[17]
2.6557
E
-
03
–
0.82
E
+
00
–
[18]‡
2.65567
E
-
03
n.a.
2.66
E
+
04
(67 iter.)
[19]†
2.6557
E
-
03
n.a.
1.11
E
+
01
(182 iter.)
♭Result with relative convergence tolerance of 1.0E-02. Computer HP J2240. §PC: Pentium IV, 2.4 GHz, 256 MB RAM. ‡Ultra SPARC-II CPU (2×360 MHz), 512 MB RAM. Optimality margin ϵ=1E-02. †Relative convergence tolerance of 1.0E-03. PC: Intel Pentium 4 with 3.2 GHz.Table 13
Comparative results forIRR1 (ABC/200n).
Source
J
b
e
s
t
J
a
v
g
T
a
v
g
N
F
E
a
v
g
1.185851
E
-
06
2.381824
E
-
03
1.21
E
+
01
449
[4]
1.18584
E
-
06
n.a.
–
–
[10]♭
1.1858
E
-
06
n.a.
1.08
E
+
01
(38 iter.)
[16]§
1.185845
E
-
06
1.185845
E
-
06
0.46
E
+
00
–
[17]
1.1858
E
-
06
–
0.33
E
+
00
–
[18]‡
1.18562
E
-
06
n.a.
7.67
E
+
02
(37 iter.)
[19]†
1.1858
E
-
06
n.a.
0.23
E
-
01
(4 iter.)
♭Result with relative convergence tolerance of 1.0E-02. Computer HP J2240. §PC: Pentium IV, 2.4 GHz, 256 MB RAM. ‡Ultra SPARC-II CPU (2×360 MHz), 512 MB RAM. Optimality margin ϵ=1E-02. †Relative convergence tolerance of 1.0E-03. PC: Intel Pentium 4 with 3.2 GHz.Table 14
Comparative results forIRR2 (ABC/200n).
Source
J
b
e
s
t
J
a
v
g
T
a
v
g
N
F
E
a
v
g
4.411030
E
-
06
4.593098
E
-
06
1.55
E
+
01
447
[26]†
1.7604
E
-
04
1.7604174
E
-
04
–
–
†Best solution obtained after 100 multistarts running Dynopt.Table 15
Comparative results forREV_A (FA/200n).
Source
J
b
e
s
t
J
a
v
g
T
a
v
g
N
F
E
a
v
g
1.933663
E
-
07
2.119351
E
-
06
3.17
E
+
01
953
[10]†
3.367
E
-
07
n.a.
2.73
E
+
02
(56 iter.)
†Result with absolute convergence tolerance of 1.0E-08. Computer HP J2240.Table 16
Comparative results forREV_B (FA/50n).
Source
J
b
e
s
t
J
a
v
g
T
a
v
g
N
F
E
a
v
g
1.587466
E
-
03
1.617494
E
-
03
1.27
E
+
01
374
[10]†
1.586
E
-
03
n.a.
5.68
E
+
02
(349 iter.)
[19]‡
1.5875
E
-
03
n.a.
2.62
E
+
03
(40552 iter.)
†Result with relative convergence tolerance of 1.0E-02. Computer HP J2240. ‡Relative convergence tolerance of 1.0E-03; PC/Intel Pentium 4 with 3.2 GHz.Table 17
Comparative results forTHEO1 (HS/50n).
Source
J
b
e
s
t
J
a
v
g
T
a
v
g
N
F
E
a
v
g
1.070635
E
-
03
1.090058
E
-
03
9.11
E
+
00
255
[26]†
4.01
E
-
03
4.0303
E
-
03
–
–
†Best local solution after 100 multistarts running Dynopt.Table 18
Comparative results forTHEO2 (ABC/200n).
Source
J
b
e
s
t
J
a
v
g
T
a
v
g
N
F
E
a
v
g
3.578108
E
-
02
4.220937
E
-
02
2.12
E
+
01
448
[26]†
2.595
E
-
02
2.595
E
-
02
–
–
†Solution obtained after 100 multistarts running Dynopt.Table 19
Comparative results forVOL_B (DE/200n).
Source
J
b
e
s
t
J
a
v
g
T
a
v
g
N
F
E
a
v
g
8.667010
E
-
01
9.384865
E
-
01
7.89
E
+
01
440
[10]†
1.3194
E
-
03
n.a.
1.00
E
+
05
(1000 iter.)
[19]‡
1.2492
E
-
03
n.a.
4.30
E
+
01
(536 iter.)
†Computer HP J2240. ‡Relative convergence tolerance of 1.0E-03; PC/Intel Pentium 4 with 3.2 GHz.Overall, after all the numerical comparisons, it is possible to conclude that the selected metaheuristics when enhanced with the local intensification phase are effective in achieving good quality solutions with a reduced computational effort. Furthermore, we have also shown that the sequential direct method combined with the selected metaheuristics competes very favorably with exact methods and other metaheuristics available in the literature.


## 5. Conclusions

In this paper, we have analyzed the performance of five well-known metaheuristics when solving parameter estimation problems in dynamic system models. The sequential numerical direct method is applied to the DMbPE problem and the resulting optimization is performed directly making use of a numerical integration formula for solving the system of ODE. Using nine DMbPE problems and different experimental data, error-free and random error-added data, a total of 12 instances have been used in the comparative experiments. The solutions produced by the metaheuristics have been analyzed in terms of quality, by stopping the algorithms after a specified number of function evaluations, and compared using statistical hypotheses testing. The statistical tests show that the average obtained solutions and the average CPU time are considered to be mostly significantly different, at a significance level of 5%. The best solutions produced by the metaheuristics are analyzed by means of the performance profiles. After the graphical comparisons for different numbers of allowed function evaluations, we are able to conclude that the DE is the metaheuristic that has the highest probability of giving the lowest value of the objective functionJ mostly. It is followed by HS. The FA, ABC, and DTS metaheuristics are, in this sequence, the least effective. Thus, our recommendation for a reader that favors good quality solutions falls in the DE and HS metaheuristics. On the other hand, ABC and DE are in general the metaheuristics that require the least computational effort, meaning that, for a limited computational budget, our recommendation is to choose one of them. Finally, from the comparisons with other results in the literature, involving both sequential and simultaneous direct methods and stochastic as well as deterministic global optimization methods, we state that the sequential direct method combined with one of the above recommended metaheuristics is a good alternative to solve DMbPE problems. Good quality solutions with a reduced computational effort are thus provided.



---

*Source: 3213484-2018-01-29.xml*