Real-time magnetic dipole detection with single particle optimization

Applied Soft Computing 23 (2014) 460–473

Contents lists available at ScienceDirect

Applied Soft Computing journal homepage: www.elsevier.com/locate/asoc

Real-time magnetic dipole detection with single particle optimization Giovanni Iacca ∗ , Frank Lennart Bakker, Heinrich Wörtche INCAS3 , P.O. Box 797, 9400 AT Assen, The Netherlands

a r t i c l e

i n f o

Article history: Received 11 March 2013 Received in revised form 9 June 2014 Accepted 19 June 2014 Available online 9 July 2014 Keywords: Magnetic dipole detection Single particle optimization Online optimization Real-time optimization Computational intelligence

a b s t r a c t In recent years, the use of magnetic field measurements has become relevant in several applications ranging from non-invasive structural fault detection to tracking of micro-capsules within living organisms. Magnetic measurements are, however, affected by a high noise due to a number of causes, such as interference from external objects and the Earth magnetic field. Furthermore, in many situations the magnetic fields under analysis are time-variant, for example because generated by moving objects, power lines, antennas, etc. For these reasons, a general approach for accurate real-time magnetic dipole detection is unfeasible, but specific techniques should be devised. In this paper we explore the possibility of using multiple 3-axis magnetic field sensors to estimate the position and orientation of a magnetic dipole moving within the detection area of the sensors. We propose a real-time Computational Intelligence approach, based on an innovative single particle optimization algorithm, for solving, with an average period of 2.5 ms, the inverse problem of magnetic dipole detection. Finally, we validate the proposed approach by means of an experimental setup consisting of 3 sensors and a custom graphical application showing in real-time the estimated position and orientation of the magnetic dipole. Experimental results show that the proposed approach is superior, in terms of detection error and computational time, to several state-of-the-art real-parameter optimization algorithms. © 2014 Elsevier B.V. All rights reserved.

Introduction For many industrial applications it is desirable to obtain an accurate characterization of a ferromagnetic object (or a set of objects), by measuring the magnetic fields induced by it. Magnetic field detection can be used, for example, for identifying vehicles based on their magnetic profile, or for detecting defects, e.g. inner cracks or micro-holes, in metallic structures (rails, steel beams, etc.). Broadly speaking, all those applications which rely on the analysis of magnetic field measurements require that the magnetic dipoles (or in general, magnetic domains) present in the structure under study are precisely localized. In addition to that, in some cases it might be needed to detect a time-variant magnetic field in real-time, e.g. when one wants to track a magnetic dipole whose position and orientation change with time. A paradigmatic example of a real-time application of this kind would be a magnetic-based motion tracking system, where one or more dipoles are attached to a moving object and their position/orientation must be estimated online. Such a system could be used, for instance, in video gaming, motion capture, 3D painting, or interactive path-planning of industrial robots.

∗ Corresponding author. Tel.: +39 3898088796. E-mail addresses: [email protected], [email protected] (G. Iacca), [email protected] (F.L. Bakker), [email protected] (H. Wörtche). http://dx.doi.org/10.1016/j.asoc.2014.06.026 1568-4946/© 2014 Elsevier B.V. All rights reserved.

Despite the broad potential applications, the problem of detecting and interpreting magnetic fields in real-time still represents an open challenge. The reasons are manifold, namely: (1) the high noise which usually affects magnetic measurements; (2) the effect of the static earth magnetic field; (3) the interference due to the magnetic fields generated by surrounding electrical currents and objects, such as computers and other electronic devices; (4) the real-time constraints imposed by the sensor data acquisition frequency and, in case of moving magnetic dipoles, their speed. In recent years, several methods have been developed for localizing magnetic dipoles, and have been applied to various contexts ranging from engineering to biomedicine. These methods are based, in general, on the solution of the inverse magnetic problem by means of different optimization strategies, such as Computational Intelligence algorithms or classic local search methods. In [1], an array of 3-axis magnetic sensors is used to track a magnet, modeled as a single dipole, inside a capsule endoscope for the interactive diagnosis of the gastrointestinal tract; in this case, the localization is solved by means of a combination of a linear localization algorithm [2] and the Levenberg–Marquart method [3], where the first provides a valid initial guess and the latter tries to improve upon it. In [4], the same authors extend this approach for detecting position and orientation of a rectangular magnet in a capsule endoscope by means of basic genetic algorithm (GA) [5] and particle swarm optimization (PSO) [6]. In both studies it was obtained, in

G. Iacca et al. / Applied Soft Computing 23 (2014) 460–473

simulation, a high localization/orientation accuracy and computational time in the order of tenths of a second – up to 2 s – depending on the desired accuracy. Paper [7] proposes a similar approach for tracking, by means of an adhesive magnet, the tongue movements within the 3-D oral space; in this case PSO was compared with three local search methods, namely the DIRECT [8], Powell [9], and Nelder–Mead [10] algorithms, with the Powell algorithm resulting in the highest accuracy and an average processing time of 43.9 ms/sample. Some alternative approaches, which do not make use of an explicit optimization algorithm, have also been proposed in the literature. For example, in [11], a real-time localization technique is described, based on the Euler homogeneity equation. In [12], a Hopfield neural network [13] was applied to the estimation of magnetic moments, while a self-adaptive evolutionary algorithm was used for training the model by minimizing the least-squares error between the observations and the data produced by the neural network. Also in these cases, the proposed approaches display a high accuracy, a good robustness, and a relatively simple implementation. Most of the methods from the literature, however, either make problem-specific assumptions or they are not suitable for fast realtime tracking of moving magnetic objects (i.e. in the order of few milliseconds per sample), because for instance they require an offline learning phase or an offline data processing. To fill this gap, in this paper we present a novel optimization-based method for fast detection of position and orientation of a single magnetic moment in space. The proposed optimization algorithm, named SPORE, is characterized by a single particle which explores the parameter space according to an extremely simple self-adapting search logic. Thanks to its simplicity, SPORE is way faster than other methods previously proposed in the literature, being able to obtain sufficiently good results in a remarkably short computational time (on average, 2 ms, that is from two to three orders of magnitude smaller than other methods in the literature). Additionally, since our method employs a single particle at every step of the optimization process, its memory footprint is dramatically lower than, for example, population-based algorithms such as GA or PSO (as used in [4]), or other modern algorithmic structures employing covariance matrices or solution archives. Thus our method not only is suitable for fast real-time magnetic moment detection, but also allows for possible implementations on embedded systems endowed with limited CPU power and memory. This might be, for example, the case of portable biomedical detectors or portable magnetic scanners that could be used in structural non-invasive inspection. To test the proposed approach, we built an experimental setup composed of three off-the-shelf 3-axis magnetic field sensors, and designed a Java Graphical User Interface (GUI) in which the position and orientation of the magnetic moment is displayed, in a 3-D environment, on the screen of a computer connected to the sensors via USB. An optimization process is performed in background every time a new set of sensor data (3 samples X–Y–Z, one per sensor) is available for fitting, and the GUI is updated in real-time while the magnetic moment moves around the sensors’ detectable area. To asses the superiority of our method in terms of computational time and fitting error, we compared it with a selection of state-ofthe-art real-parameter optimization algorithms, including recent variants of PSO and differential evolution (DE), as well as some modern single solution optimization algorithms and classic local search methods. The remainder of this paper is organized as follows. The next section describes the physical model of a single dipole moment and how its inverse formulation can be interpreted as a real-parameter optimization problem. Section ‘Experimental setup’ describes the experimental setup used for tests, while in Section ‘Single Particle

461

z

θ

m

y φ

x Fig. 1. The orientation of m is represented in spherical coordinates (M, , ).

Optimization with Restarted Exploration’ we describe our new single particle optimization algorithm. Section ‘Experimental Results’ presents the experimental results obtained with SPORE in comparison with other real-parameter optimization algorithms. Finally, Section ‘Conclusions and future works’ gives the conclusion of this study and suggests a few possible industrial applications that might be investigated in the future. Background: magnetic dipole moment In order to find the location, orientation and spatial structure of magnetic objects simply by measuring the magnetic fields they induce, one needs to be able to find the solution to an inverse problem. Although the magnetic field around a magnetic object can be calculated relatively easily, doing the inverse operation (i.e. finding its structure, position and orientation based on the magnetic field measurements) is generally more difficult, as it requires the use of an efficient fitting algorithm and often also a list of assumptions to limit the search space. Similarly to what was done for instance in [1] and [2], in this section we first formulate the direct expression of the external magnetic field induced by a single magnetic dipole moment, and then we define its inverse problem as an optimization problem. Direct formulation The external magnetic field H(r) induced by a single magnetic moment m is given by [14]: H(r) =



1 4||r||3

3r(m · r) ||r||2



−m

(1)

where r denotes the vector from the magnetic moment to the point in space where the field is calculated. Because the magnitude of the magnetic moment is usually fixed and for symmetry reasons, we choose for a representation of the moment orientation in spherical coordinates. Fig. 1 shows the general convention we use for the coordinate system of m. The angle  represents the angle between the Z-axis and the XY-plane, whereas  denotes the angle within the XY-plane. Hence, the orientation of m in Cartesian coordinates can be rewritten as:

⎛

M sin() cos()

m(x, y, z) = ⎝ M sin() sin()

⎞ ⎠

(2)

M cos() where M = ||m|| is constant for a single dipole moment. The position of the magnetic moment is denoted by the vector rm (x, y, z). The magnetic field present at the position of the sensors is then given by Eq. (1) for r = rs − rm where rs (x, y, z) is the position of the magnetic field sensor.

462

G. Iacca et al. / Applied Soft Computing 23 (2014) 460–473

Inverse problem When the sensor positions rs (x, y, z) are known precisely, and combining this information with the magnetic fields measured at the sensor positions, the inverse problem associated to the direct problem in Eq. (1) can be solved by using a fitting method. Indicating with p the parameter vector [rm (x, y, z), M, , ] describing the unknown magnetic dipole, the magnetic data fitting process (i.e. the magnetic moment detection problem) can be formulated as follows: find

p∗

s.t.

f (p∗ ) = min f (p) p ∈ D

(3)

where D is the search space and * indicates the optimal solution in D. With reference to Eq. (1), the objective function f(·) (in the following, simply referred to as “fitness”1 ) is a measure of the fitting error, that can be expressed for instance in the form of sum of squared errors:

 

Nsensors 3

f (p) =

i=1

˜ j (p, rs (xi , yi , zi ))]2 [Hi,j − H

(4)

j=1

where i indicates the sensor index and j the Cartesian component index (j = 1, 2, 3 for X, Y and Z respectively); Hi,j represents the jth Cartesian component of the magnetic field measured at the ˜ j (p, rs (xi , yi , zi )) is the jth Cartesian component of the ith sensor; H (estimated) magnetic field generated by a magnetic moment characterized by the parameters p = [rm (x, y, z), M, , ] which would be detected at the sensor position rs (xi , yi , zi ). As said these positions rs (xi , yi , zi ), i = 1 . . . Nsensors , must be known a priori before starting the optimization process and are considered constant. It can be easily seen that lower values of f(·) indicate that the magnetic moment is detected correctly. Theoretically, a perfect detection would produce a fitting error equal to zero; in practice, considering the noise on the sensor data and some numerical errors due to the sensor calibration process (see below), it is impossible to obtain such a value but rather one should minimize the error to achieve a satisfactory fitting. As a final remark, we should note that with three magnetic sensors the inverse problem defined in Eq. (3) has a unique solution, since with 9 measurements (X–Y–Z per each 3-axis sensor) and 6 unknown parameters [rm (x, y, z), M, , ], this problem is equivalent to solving an over-determined non-linear system of 9 equations in 6 variables. This is also true for any number of sensors larger than three. On the other hand, with only one sensor we would have an infinite number of solutions (6 unknowns, 3 measurements). With two sensors, we would have a limit situation of a non-linear system with 6 unknowns and 6 measurements, which might not be univocally solvable due to the spacial symmetry of some magnetic field configurations. We should point out that this analysis is valid only if we are interested in determining position, orientation and magnetic moment M. With a lower number of unknowns, the above analysis differs. Experimental setup In order to test and verify the idea of the single dipole localization, an experimental demonstrator setup has been realized. The setup consists of three 3-axis magnetic RM3000 field sensors [15]. These sensors provide a continuous flow of measures of the

1 We should remark that this terminology is inherited from classic evolutionary and swarm intelligence algorithms, where each solution of a problem is seen as an “individual” and its performance with respect to the optimization problem is seen as its “fitness”.

Fig. 2. Setup with three 3-axis magnetic sensors. The sensor absolute positions refer to a coordinate reference system whose origin is the center of the figure. Axis lines X and Y are oriented as shown by the arrows, while the Z-axis points up. For the demonstrator setup, a magnetic moment will be detected within a box of 40 cm × 40 cm × 20 cm (see Section ‘Experimental results’).

magnetic field strength in the X, Y and Z direction. The sensor count-rate was set to 500 counts in standard mode, providing a gain of about 100 counts per ␮T, a sample frequency (per axis) of about 200 Hz, and a minimum noise level of 10 nT. Each sensor is mounted on a cased custom PCB developed in-house, and connected via an FTDI cable to an ATmega8U2 Breakout board which acts as FTDI/USB bridge. The three breakout boards are connected to a USB mini-hub which sends the sensor data to the PC running the real-time demonstrator described in Section ‘Real-time demonstrator’. The top-view of the experimental setup is shown in Fig. 2. The three sensors are mounted on an aluminum plate and the sensor positions are measured with an accuracy of about 5 mm. We should remark that each sensor actually contains 3 independent sensing elements, one for each axis, which are mounted at different positions on the sensor PCB. For the optimization procedure, however, we assume only a single position rs (xi , yi , zi ), see Eq. (4), located on the center of symmetry of the PCB case. This introduces a small error in the localization of the magnetic moment. Another source of error is the custom PCB board/FTDI interface, which due to the absence of a power ripple filter introduces an additional background noise in the order of few tens of nT. Commercial high-quality boards should be used in industrial contexts to reduce this noise component.

Sensor calibration In order to use the sensor data for fitting, it is important to ensure that those data are calibrated, i.e. their order of magnitude and orientation are consistent with the real magnetic field induced by the dipole. In our experiments, a careful calibration procedure has been carried out, before installing the sensors on the aluminum plate, using the approach described by Camps et al. [16]. The calibration

G. Iacca et al. / Applied Soft Computing 23 (2014) 460–473

procedure for a 3-axis sensor is based on the assumption that its output can be written as:

⎛

Vx (t)

⎞

⎛

˛x

⎜ ⎟ ⎜ ⎝ Vy (t) ⎠ = ⎝ sxy Vz (t)

sxz

⎞⎛

˛y

syz ⎠ ⎝ hy (t) ⎠ + ⎝ ˇy ⎠

syz

˛z

⎟ ⎜

hz (t)

ˇx

⎞

sxz

⎟⎜

hx (t)

⎞ ⎛

sxy

⎟

(5)

ˇz

where V(t) are the raw sensor measures, ˛ are the sensor gains, s are the non-orthogonality factors (which take into account the misalignment between the coordinate reference system of the sensors and the absolute reference system of the real magnetic field), h(t) are the real magnetic field and ˇ are the sensor offsets. To calibrate each sensor, we collected the sensor readings while slowly rotating the sensor (in every possible direction) in the Earth magnetic field, i.e. in absence of any other magnetic moment dipole. In total, we collected 950 different sensor values for each sensor. By assuming:



||H|| =

hx (t)2 + hy (t)2 + hz (t)2 = K

(6)

where K is the constant, known, Earth magnetic field (to compute its value at our location, 53◦ 0 N 6◦ 33 E, we used the International Geomagnetic Reference Field 2010 model, see [17]) one can fit the recorded data with Eq. (5) to find the unknown parameters ˛, s and ˇ. As in [16], we used the Levenberg–Marquardt algorithm to minimize a quadratic fitting error function. For each sensor, we performed this procedure one-time-only and then saved the calibration parameters into a configuration file. Such file is read before starting the online optimization process (see below), so that the raw measures V(t) coming from the sensors are converted into the corresponding calibrated measures h(t) by applying Eq. (5). Online optimization process Algorithm 1.

Online optimization process

// initialization initialize cyclescal read the calibration parameters ˛, s and ˇ from the configuration file open and configure sensor USB communication // measure static field for i = 1 → cyclescal do HSraw ← readFromSensors() end for ¯ Sraw compute average static field H ¯ S ← calibrate H ¯ Sraw according to Eq. (5) H cal

// online optimization initialize ε, threshold, Dlarge , Dsmall , Tlarge , Tsmall while user does not stop do Hraw ← readFromSensors() Hcal ← calibrate Hraw according to Eq. (5) ¯S H ← Hcal − H cal

if p* = / NULL then algorithm.initialize(Dsmall , Tsmall , p* ) else algorithm.initialize(Dlarge , Tlarge ) end if p* , f(p* = algorithm.run(H, rs (xi , yi , zi ), i = 1 → Nsensors ) if f(p* ) < ε then if M > threshold then update GUI (update dipole position/orientation) and log results else update GUI (remove dipole) p* = NULL end if end if end while

An overview of the online optimization process is given in Algorithm 1. At the beginning of the procedure, the calibration parameters ˛, s and ˇ are read from the configuration file (as described in the previous section) and the USB channel is established and configured to obtain the data from the sensors. Then,

463

the static Earth magnetic field (i.e. the field measured in absence of any magnetic moment) HSraw is measured and averaged over a fixed number of cycles cyclescal . Here we indicate with readFromSensors() the routine that reads the raw data coming from the ¯S sensors. The average static field H raw is then calibrated according ¯ S , compensated for sensor gains/offsets and to Eq. (5) to obtain H cal non-orthogonality of the reference systems. In particular, to obtain the calibrated data, the linear system Eq. (5) must be solved, with h(t) unknown and V(t) equal to the raw sensor data. At the end of the initialization stage, the online process is started. This process is continued until the user stops it from the graphical interface (see next subsection). At each step of the loop, the raw sensor data HSraw are collected and calibrated as described previously, to obtain HScal . In addition to that, the average static (calibrated) ¯ S is subtracted from the calibrated data HS to obtain Earth field H cal

cal

the magnetic field H that has to be fitted. For the sake of clarity, we should note that all the H variables introduced so far are implemented as matrices in RNsensors ×3 , where Nsensors is the number of sensors (in our case, 3). Once H is computed, an optimization algorithm (either the proposed SPORE or one of the comparison methods, see Sections ‘Single particle optimization with restarted exploration’ and ‘Experimental results’) is executed in the attempt of solving the optimization problem in Eq. (3), i.e. fitting H and detect Cartesian position, magnitude and orientation of the unknown magnetic moment (p* ). With reference to Algorithm 1, the functions algorithm.initialize() and algorithm.run()2 correspond, respectively, to the parameter initialization and the execution of the optimization algorithm. If the fitting error f(p* ) is smaller than a given threshold ε and the detected moment is sufficiently large (i.e., bigger than a given threshold), the fitting process is considered successful and its output is provided (either graphically or on a log file) to the user. It should be noted that, in order to provide a plausible initial guess for each step of the algorithm, at the beginning of each iteration the algorithm is initialized with the best solution p* found in the previous step. In case of population-based optimization algorithms, this is equivalent to the injection of p* in the initial randomly sampled population; in case of single solution algorithms, this means that the initial solution is just p* . If the detected moment is smaller than the threshold, however, the previous p* is nullified, since this means that no magnetic moment is present in the detectable area and that the previous solution shouldn’t be retained anymore. When p* is null, the algorithm explores the whole search space Dlarge (spanning over a large detectable area above the aluminum plate) using a larger computational budget, Tlarge ; otherwise, the search space is shrunk over a smaller portion of the search space, Dsmall , centered around the previous best solution p* . In this case the algorithm is also allotted a shorter budget, Tsmall , under the assumption that a shorter number of fitness evaluations is needed to converge towards the optimum when a “good” initial guess is provided. This mechanism is especially useful if the magnetic moment moves relatively smoothly into the search space, so that only incremental variations of its position and orientation must be detected. On the other hand, when no previous initial solution p* is available, the algorithm must be allowed to search the whole detectable area, so to estimate any new dipole position/orientation within it. Real-time demonstrator In order to display the estimated position and orientation of the dipole moment, we developed a multi-platform Java SWT

2 We use here an object-oriented notation to highlight the fact that these methods refer to each single optimization algorithm.

464

G. Iacca et al. / Applied Soft Computing 23 (2014) 460–473

Fig. 3. Graphical user interface of the demonstrator setup. The interface shows the position (dot) and orientation (line) of the magnetic dipole moment and includes the possibility of tracking the moment in real-time. The three dots on the XY-plane correspond to the positions of the three sensors shown in Fig. 2.

application based on the OpenGL library JOGL [18,19]. The USB connection is handled using the JSSC library [20], while JAMA [21] is used for solving the calibration linear system in Eq. (5). The software runs, in background, the online optimization process shown in Algorithm 1 and updates, in real-time, a 3-D graphical environment showing the estimated dipole moment. Two update modes are available, namely track, which tracks the estimated position/orientation trajectory, and on-the-fly, which only plots the estimated position and orientation obtained from the last sensor samples. Through the interface, the user can control the optimization process (pausing or stopping it), set the update mode, and enable the visualization of the estimated position/orientation of the dipole and the sensor positions. Fig. 3 shows a freehand trajectory of a magnetic dipole moment where both the position and orientation of the moment are determined using the procedure illustrated in Algorithm 1 with the single particle optimization algorithm described in the following section. Single particle optimization with restarted exploration In the last two decades, many researchers have put a great effort in designing robust and reliable optimization algorithms capable of guaranteeing good results on broad sets of benchmark functions and real-world problems. Nowadays, several state-of-the-art optimization methods can be rightfully considered general-purpose optimizers. Examples in the area of Evolutionary Computation and Swarm Intelligence are represented by CMA-ES [22] and various modifications of it [23–26], as well as numerous variants of Differential Evolution [27–35], Particle Swarm Optimization [36–44], and Memetic Computing approaches [45–52]. However, despite the availability of these powerful optimization methods, there are still some very specific real-world optimization problems (see e.g. the robotic and sensor applications presented in [53–60]) which are difficult to tackle by means of general-purpose optimizers. For these cases, it is still reasonable to devise ad hoc optimization methods, or adapt existent algorithms to the problem, for instance embedding domain-specific knowledge into the algorithmic operators.

In our case study, the real-time constraints and the presence of noise on the sensor measurements pose some technical challenges which make the optimization process particularly hard. In order to perform real-time detection of a magnetic dipole, the optimization algorithm should have the following features: (a) low computational overhead, i.e. it should perform a minimal number of mathematical operations; (b) strong exploitation, i.e. it should be able to obtain good results in a very limited number of fitness evaluations; (c) ability to escape local optima, since the fitness landscape corresponding to Eq. (4) is characterized by multiple “peaks”, mainly due to the noise and to the non-linear sensitivity of sensors. Having in mind these goals, we propose here a novel, fast single particle optimization named SPORE (Single Particle Optimization with Restarted Exploration). In a nutshell, SPORE consists of a single particle that is repeatedly perturbed, variable-wise, as long as the perturbations produce a fitness improvement. This mechanism endows the algorithm with a strong exploitation pressure. If the perturbations are unsuccessful, a simple restart scheme is used to partially re-sample the particle, inheriting a variable-size portion of the current best solution. The latter component gives the algorithm the ability of escaping local optima, still preserving part of the best solution found so far. This process goes on until a stop condition is met. Let us analyze more in detail the working principles of SPORE. The pseudo-code of the algorithm is shown in Algorithm 2, where it can be seen that the method consists of an initialization phase, followed by an optimization loop. Each cycle of the loop consists in turn of two steps, namely: (1) a variable-wise perturbations loop and (2) a restart with partial inheritance. During the initialization, a parameter ıi is initialized for each variable to the width of the search space along that dimension (ubi − lbi , where ubi and lbi are, respectively, the upper and lower bound of the search space along the ith direction). Together with that, for each variable a maximum of perturbation iter number

ations is initialized, kmax = log2 ıi / , where  is an algorithm i

parameter which controls the minimum perturbation range. After the initialization, the optimization loop starts. At the beginning of the first step of the loop, the current best fitness function is saved. Then, each ith variable is perturbed repeatedly, starting with a perturbation step ı = ıi /2. More specifically, for each ith variable, and for a maximum number of perturbation cycles (kmax ), the following operations are performed: first, a perturbai

tion of radius ı is attempted in one direction; if this move leads to an improvement, the radius ı is halved and another perturbation in the same direction is attempted, otherwise the ith variable is reset to its previous value and the opposite direction −ı is used in the next perturbation cycle. When for the current ı both directions have been tried and no further improvements are possible, ı is halved and a new perturbation cycle begins. This way incrementally smaller perturbations are attempted, until the maximum number of perturbation cycles kmax is reached. In order to introi

duce the possibility of interrupting this exploitative search when the best solution is not improved, at the end of each cycle the per2 turbation loop is stopped with probability (k/kmax ) , where k is i

the index of the cycle. In this way, the perturbation loop has a very low chance of being interrupted during the first cycles, while this chance increases with the number of unsuccessful perturbation cycles. Whenever the first step does not produce any fitness improvement, the second step, i.e. restart with partial inheritance, is activated. This mechanism is introduced in the algorithmic structure in order to prevent the algorithm from being trapped into local optima, and balance the exploitation pressure of the first step. As shown in Algorithm 3, a logics similar to the binomial crossover

G. Iacca et al. / Applied Soft Computing 23 (2014) 460–473

used in differential evolution [61] is adopted. At restart activation, a dynamically adapting crossover rate CR is computed based on the current and maximum number of fitness evaluations, respectively neval and Neval . The crossover rate is updated in such a way that it takes its maximum value at the beginning of the optimization, and decreases during the next stages. Once CR has been computed, for each ith variable of the current best particle a uniform random number r in [0, 1] is generated and compared with it. If r < CR, the ith variable is re-sampled within its bounds [lbi , ubi ]; otherwise it will take the value of the corresponding variable in the current best solution. The modified particle x is then evaluated, and, in case of fitness improvement, it is copied into the current best particle. The main idea behind SPORE is essentially that, in order to address the optimization problem at hand, an extremely exploitative component is needed to refine the search within the current basin of attraction. Similar to the Hooke–Jeeves pattern search [62], the variable-wise perturbations loop tries to improve upon the current solution applying incrementally smaller perturbations, in both directions, on each variable. A randomized stop criterion is introduced within this component so to allow a higher chance of escaping local optima. The rationale for the randomized stop is that the algorithm should perform a larger number of perturbation cycles when exploiting a very promising region of the search space (thus obtaining, with a high probability, incremental improvements), and a small number of perturbations when it is stuck into a local optimum, where it fails at finding further improvements. Therefore, introducing a randomized number of perturbations (and a stop probability which quadratically increases with it), we endow the algorithm with a very simple form of adaptation, and, therefore, with a greater flexibility. It should be noted that the quadratic function guarantees a smoother probability increase than a linear 2 function (because (k/kmax ) ≤ (k/kmax ) ∀k ∈ [1, kmax ]). i

i

i

On top of that, a simple restart mechanism is used to resample the current particle whenever the exploitative component fails at improving upon it. This re-sampling mechanism guarantees a better trade-off between exploration and exploitation in two ways: on one hand, it retains parts of the current particle (rather than sampling a solution ex novo), thus preserving part of the previous achievements; on the other, the adaptive crossover rate adjusts the amount of variables which are inherited from the best solution, so that at the beginning of the optimization only a small number of variables are inherited (stronger exploration pressure), while subsequently larger parts of the best solution are retained (thus gradually favoring the exploitation pressure). Algorithm 2.

SPORE pseudo-code

// initialization initialize initial particle x initialize best best = x fx = fbest = f(x) initialize  for i = 1 → n do initialize ıi = (ub i − lbi ) kmax = log2 ıi / i

end for initialize direction flag dir = 0 // optimization loop while stop criterion is not met do // step 1: variable-wise perturbations loop save current best fitness finit = fbest for i = 1 → n do update perturbation ı = ıi /2 stop = FALSE for k = 1 → kmax and not stop do i

save old position xbak = xi save old fitness fbak = fx update position xi = xi + ı

465

update direction flag dir = dir + 1 calculate fitness fx = f(x) if fbak < fx then change perturbation direction ı = − ı reset old position xi = xbak reset old fitness fx = fbak else halve perturbation ı = ı/2 reset direction flag dir = 0 k=k+1 end if // if both directions have been tried if dir = =2 then halve perturbation ı = ı/2 reset direction flag dir = 0 k=k+1 end if if fx < fbest then best = x fbest = fx else



if rand(0, 1) <

2

k/kmax

then

i

stop = TRUE end if end if end for end for // step 2: restart with partial inheritance if (finit − fbest ) = =0 then perform restart as in Algorithm 3 end if end while

Algorithm 3.

Restart with partial inheritance pseudo-code

adaptive crossover rate CR = [(Neval − neval )/Neval ]2 for i = 1 → n do if rand(0, 1) < CR then xi = rand(lbi , ubi ) else xi = besti end if end for calculate fitness fx = f(x) if fx < fbest then best = x fbest = fx end if

Experimental results Several optimization algorithms have been considered, besides SPORE, in order to minimize the fitting error in Eq. (4). During our experimental campaign, we tested both modern global optimization algorithms and classic local search methods, namely: • Modern population-based algorithms – Self-Adaptive Differential Evolution + pBX crossover (MDEpBX) [63], with population size equal to 100 individuals and group size q equal to 15% of the population size. – Comprehensive Learning Particle Swarm Optimizer (CLPSO) [37], with population size equal to 60 individuals. – Covariance Matrix Adaptation Evolution Strategy (CMA-ES) [22], with the default parameter setting of the original Java implementation [64], namely  = 4 + 3 ln(n) ,  = /2 , and initial step-size = 0.2. • Modern single solution algorithms – 3 Stage Optimal Memetic Exploration (3SOME) [65], with inheritance factor ˛e = 0.05, middle distance exploration hypercube size ı equal to 20% of the total decision space width, coefficient of generated points at each activation of the

466

G. Iacca et al. / Applied Soft Computing 23 (2014) 460–473

middle distance exploration k = 4, short distance exploration radius = 0.4 and local budget fixed to 150 iterations. – Intelligence Single Particle Optimization (ISPO) [66], with acceleration A = 1, acceleration power factor P = 10, learning coefficient B = 2, learning factor reduction ratio S = 4, minimum learning threshold E = 1e−5, and learning steps PartLoop = 30. – non-uniform Simulated Annealing (nuSA) [67], with mutation exponent B = 5, temperature reduction ratio ˛ = 0.9, temperature reduction period Lk = 3, and number of initial solutions to set the initial temperature initialSolutions = 10. • Classic local-search algorithms – Powell Algorithm [9], using Brent method [68] as (1dimensional) line optimizer, with 100 fitness evaluations per each line. – Nelder–Mead Algorithm [10], with reflection coefficient ˛ = 1, contraction coefficient ˇ = 0.5, expansion coefficient  = 2 and shrinkage coefficient ı = 0.5. – Rosenbrock Algorithm [69] with positive perturbation factor ˛ = 2, negative perturbation factor ˇ = −0.5, and threshold for coordinate system rotation = 10−5 . – the quasi-Newton algorithm implemented in the package UNCMIN [70], with finite difference approximation of gradient and Hessian, and default parametrization apart from tolerance on gradient and step size set to zero. For each algorithm, we used the parameter setting suggested in the original paper. As for SPORE, we set the parameter  to 10−8 . To handle the search space bounds, we implement in all the algorithms a toroidal mechanism, consisting of the following: given an interval [a, b], if xi = b + , i.e. the ith design variable exceeds the upper bound by a quantity , its value is replaced with a + . A similar mechanism is applied for the lower bound. As shown in Algorithm 1, we defined a “large” computational budget Tlarge = 250 × n fitness evaluations, and a “small” one, Tsmall = 60 × n, where n is the number of problem variables (6 in our case, [rm (x, y, z), M, , ]). As for the other parameters required, we set cyclescal = 100 (calibration cycles), ε = 105 (fitting error threshold), and threshold = 20 ␮T. These values were chosen empirically by performing some prior experiments. Finally, the “large” search space Dlarge was set to: [[−0.2, 0.2]; [−0.2, 0.2]; [0, 0.2]; [0, 200], [0, ], [0, 2]] where each ith pair corresponds to the lower and upper bound of the ith variable. Units are m for x, y and z, ␮T for M, and rad for  and . The “small” search space Dsmall was defined as a 6-dimensional hyper-volume, centered around the current p* , and whose side along each dimension is 1/10 of the width of the corresponding large search space. This value was chosen based on the worst-case optimization time so to bound the maximum allowed/expected speed of the dipole. We must remark that the large search space corresponds to the entire detectable area of the sensors, while the small one has been chosen under the assumption that the dipole moves smoothly and with a relatively limited speed. Simulation results In order to assess the precision of SPORE, we performed two different simulated experiments where the position/orientation trajectory of the magnetic dipole moment was generated numerically. In the first experiment, we kept the dipole orientation fixed and made the dipole follow a circular trajectory. In the second, we kept the position fixed and changed the orientation according to a periodic law. In both cases we fixed the dipole magnitude to a constant. Nevertheless, in both the simulations we used the aforementioned real-time optimization process to estimate all six

parameters [rm (x, y, z), M, , ]. The search space was defined as described above. The two trajectories can be detailed as follows: • First experiment

⎧ x ⎪ ⎪ ⎪ ⎪ y ⎪ ⎪ ⎪ ⎪ ⎨z

M

⎪ ⎪  ⎪ ⎪ ⎪ ⎪ ⎪  ⎪ ⎩

=

0.1 · sin(2f (k − d))

=

0.1 · cos(2f (k − d))

=

0.1 + 0.02 · cos(2f (k − d))

=

50

=

/2

=

0

(7)

where f = 0.002, d = 500 and k is the sample index. • Second experiment

⎧ x ⎪ ⎪ ⎪ ⎪ y ⎪ ⎪ ⎪ ⎪ ⎨z

M

⎪ ⎪  ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩

=

0

=

0

=

0.1

=

50

=

 · cos(2f (k − d))

=

 +  · cos(2f (k − d))

(8)

where f = 0.002, d = 500 and k is the sample index. A graphical representation of the two trajectories is illustrated in Fig. 4. For both the experiments, we ran 24 independent simulations for each of the eleven algorithms listed before. Each simulation was continued for 550 samples (in this case there is no calibration). Numerical results are shown in Tables 1 and 2 , respectively for the first and second experiment, where we report, averaged over 24 simulations and 550 samples, the average optimization time per sample, the average localization error Error loc = 

2 2 2 (x − x˜ ) + (y − y˜ ) + (z − z˜ ) , the average absolute error on mag˜ the average absolute error on , Error  = nitude, Error M = |M − M|, ˜ and the average absolute error on , Error  = | − |. ˜ The | − | symbol ˜ here indicates the estimated parameter. Numerical results show that in the first experiment the SPORE algorithm is the second fastest among all the algorithms under study (after the Powell method), while in the second experiment it results the fastest (the second one being the Powell method). On the other hand, population-based algorithms are, as expected, much slower (up to one order of magnitude, see MDE-pBX and CLPSO). This can be explained considering the number of operations performed by these algorithms, which involve recombination of multiple solutions and other time-consuming mechanisms, e.g. the covariance matrix adaptation in CMA-ES. In comparison, SPORE performs only simple mathematical operations involving the variable-wise perturbations and the partial inheritance restart. As for the localization error, SPORE is the second best in the first experiment (after Nelder-Mead) and the third best in the second one (after 3SOME and Powell). Regarding the error on magnitude and , single solution algorithms (including SPORE) and classic methods as Powell and Nelder-Mead clearly outperform, in both experiments, modern complex algorithms. A possible explanation for the different performance is that modern population-based algorithms are largely explorative in the very first stages of the optimization, while becoming more exploitative later on. In other words, these algorithms usually require a larger computational budget (up to several thousands of fitness evaluations) to converge to the optimum value: however, as discussed before, in this setup we are purposely using a very limited budget (from 60 × n to 250 × n), so to devise a fast real-time optimization process. On the other hand, it is interesting to note that the quasi-Newton

G. Iacca et al. / Applied Soft Computing 23 (2014) 460–473

467

Fig. 4. Simulated trajectory used in the first (left) and second experiment (right).

Table 1 Simulation results of the first experiment. Optimizer

T [ms]

Errorloc [m]

ErrorM [␮T]

Error [rad]

Error [rad]

SPORE MDE-pBX CLPSO CMA-ES ISPO 3SOME nuSA Powell Nelder-Mead Rosenbrock UNCMIN

2.192e+00 ± 3.850e−01 2.097e+01 ± 2.287e+00 1.031e+01 ± 2.950e−01 8.894e+00 ± 5.865e−01 3.221e+00 ± 7.371e−01 3.263e+00 ± 2.017e−01 8.279e+00 ± 2.380e−01 1.566e+00 ± 6.056e−01 2.676e+00 ± 3.963e−01 3.978e+00 ± 5.317e−01 3.338e+00 ± 2.032e−01

5.721e−03 ± 2.119e−03 6.883e−02 ± 1.019e−02 9.922e−02 ± 2.540e−03 9.247e−02 ± 9.759e−03 1.808e−02 ± 1.911e−02 7.954e−03 ± 2.324e−03 3.906e−02 ± 3.128e−03 1.255e−02 ± 3.028e−03 2.027e−03 ± 1.008e−03 8.622e−01 ± 2.813e−01 6.382e+01 ± 1.727e+00

1.588e+00 ± 3.719e−01 2.786e+01 ± 2.644e+00 3.766e+01 ± 1.421e+00 4.854e+01 ± 4.418e+00 5.446e+00 ± 7.173e+00 2.544e+00 ± 7.632e−01 1.777e+01 ± 1.196e+00 5.070e+00 ± 9.151e−01 8.337e−01 ± 3.442e−01 1.510e+01 ± 3.833e+00 6.246e+01 ± 2.494e+00

3.722e−02 ± 1.164e−02 5.089e−01 ± 6.466e−02 6.191e−01 ± 2.201e−02 5.789e−01 ± 4.768e−02 8.682e−02 ± 8.868e−02 6.567e−02 ± 1.821e−02 2.743e−01 ± 1.896e−02 8.668e−02 ± 1.488e−02 1.445e−02 ± 6.536e−03 1.958e−01 ± 5.850e−02 2.463e+00 ± 1.260e−01

3.379e−02 ± 2.049e−02 7.291e−01 ± 9.673e−02 1.094e+00 ± 3.600e−02 9.507e−01 ± 8.963e−02 9.944e−02 ± 1.163e−01 4.869e−02 ± 2.799e−02 4.334e−01 ± 3.628e−02 6.509e−02 ± 3.060e−02 9.954e−03 ± 1.019e−02 2.412e−01 ± 8.139e−02 1.563e+00 ± 4.136e−02

method shows a very poor performance (with errors up to four orders of magnitude larger than those obtained by SPORE): this can be explained considering that the use of gradient information might mislead the search in this case, especially due to the noise which makes the landscape extremely rugged. Finally, it can be seen that SPORE is, in the first experiment, the second best as for the error on  (after Nelder–Mead); in the second one, while CMA-ES is able to obtain the minimum error, SPORE displays a comparable performance to the other single solution algorithms (ISPO, 3SOME and nuSA). It is worthwhile mentioning that in the latter case classic local search methods tend to fail at detecting the time-variant optimum. For the sake of completeness, we report an example of fitness trend (i.e. the trend of the final value of f(·) reached, for each sample, at the end of the corresponding optimization process) obtained with the various algorithms during a single simulation of the first (Fig. 5) and second experiment (Fig. 6). Fig. 7 shows, instead, the localization trend error obtained in the same simulation of the first experiment. The graphic trend confirms visually

that SPORE largely outperforms modern population-based algorithms on both the experiments, while compared to classic local search methods and modern single solution algorithms it guarantees a much more stable behaviour along the two trajectories (see, for instance, the “spikes” that characterize the trends obtained with nuSA, ISPO, Nelder–Mead and Rosenbrock, corresponding to position/orientation configurations where the fitting error is relatively big). From this analysis we can conclude that for this specific application population-based methods are less suitable, in terms of computational time and localization/orientation error, than single solution algorithms and classic local search methods. Among these methods, Nelder–Mead and Powell are the most competitive, although they show the best performance only when the orientation is fixed, while they appear less able at detecting changes on the orientation ( and ). On the other hand, UNCMIN and Rosenbrock show the worst detection precision: according to our interpretation, this is due to the noise affecting the sensor measurements and, as a consequence, the fitness function. Since the presence of noise

Table 2 Simulation results of the second experiment. Optimizer

T [ms]

Errorloc [m]

ErrorM [␮T]

Error [rad]

Error [rad]

SPORE MDE-pBX CLPSO CMA-ES ISPO 3SOME nuSA Powell Nelder-Mead Rosenbrock UNCMIN

2.129e+00 ± 2.628e−01 1.178e+01 ± 1.395e+00 9.892e+00 ± 3.769e−01 7.579e+00 ± 1.018e+00 3.489e+00 ± 3.813e−01 2.808e+00 ± 1.752e−01 5.166e+00 ± 3.019e−01 2.279e+00 ± 3.263e−01 3.128e+00 ± 4.184e−01 2.763e+00 ± 3.568e−01 1.588e+01 ± 8.453e−01

2.465e−03 ± 9.529e−04 1.076e−02 ± 2.272e−03 5.978e−02 ± 1.681e−03 7.517e−02 ± 1.639e−02 3.601e−02 ± 1.184e−02 1.148e−03 ± 6.974e−04 2.766e−03 ± 1.110e−03 1.524e−03 ± 3.374e−04 2.534e−03 ± 1.182e−03 1.375e−01 ± 5.361e−01 6.009e+01 ± 7.868e−01

3.704e−01 ± 1.321e−01 2.950e+00 ± 6.291e−01 1.699e+01 ± 4.813e−01 3.579e+01 ± 8.092e+00 5.650e+00 ± 1.884e+00 1.924e−01 ± 1.284e−01 5.226e−01 ± 1.981e−01 2.956e−01 ± 5.581e−02 5.604e−01 ± 2.653e−01 1.115e+00 ± 3.375e−01 6.231e+01 ± 1.472e+00

2.952e−01 ± 5.928e−03 3.042e−01 ± 1.092e−02 4.072e−01 ± 1.344e−02 4.622e−01 ± 9.837e−02 4.132e−01 ± 5.458e−02 2.902e−01 ± 2.319e−03 2.972e−01 ± 4.636e−03 2.932e−01 ± 2.026e−03 2.982e−01 ± 7.835e−03 3.022e−01 ± 6.658e−03 3.207e+00 ± 1.567e−01

3.510e−01 ± 1.655e−02 1.453e+00 ± 5.088e−02 1.038e−01 ± 5.240e−02 5.712e−02 ± 2.211e−01 4.283e−01 ± 1.531e−01 3.504e−01 ± 1.273e−02 3.313e−01 ± 2.845e−02 1.547e+00 ± 1.787e−02 1.497e+00 ± 2.337e−02 1.502e+00 ± 6.475e−02 1.571e+00 ± 5.261e−02

468

G. Iacca et al. / Applied Soft Computing 23 (2014) 460–473

Fig. 5. Simulation results of the first experiment: fitness trend (log scale) of SPORE compared against modern population-based algorithms (top), modern single solution algorithms (center), and classic local search algorithms (bottom).

makes the fitness landscape rugged and therefore the use of gradient information ineffective, the performance of gradient-based algorithms (such as UNCMIN) or algorithms relying on approximations of the gradient (such as Rosenbrock, by means of rotated coordinate systems oriented along the gradient) is heavily penalized in this case. Instead, since SPORE does not rely on any gradient information, it displays an overall good performance in both the experiments, detecting with a high precision both changes on the position and orientation of the dipole. Real-world results To conclude this discussion and demonstrate the ability of the system to follow, with a high precision, the motion of the moment induced by the user, we report here the results we obtained with the experimental setup described in Section ‘Experimental setup’ on the freehand trajectory depicted in Fig. 3. The trajectory was generated moving a small cylindrical magnet (length 2 cm, radius 0.5 cm) within the detectable area defined by the search space Dlarge , and is made up of 300 samples characterized by different values of position and orientation. Due to the lack of a position/orientation

reference, in this case it is not possible to compute the localization, magnitude and orientation error as described in the previous section. On the other hand, we compare the fitness trends obtained with SPORE and the other ten competing algorithms, under the assumption that lower fitness values imply a good fitting and, ultimately, a higher precision in terms of estimated position and orientation. Fig. 8 shows the fitness trend obtained with the eleven algorithms under comparison. It can be seen that, among the selected algorithms, SPORE shows the most robust and stable behavior along the whole trajectory, obtaining, on average, a lower fitness value (i.e., a lower fitting error) with less spikes. On the other hand, population-based algorithms (especially CMA-ES) tend to obtain a poorer performance (as said, this might be explained with the larger initial exploration pressure), while the other single solution algorithms and the classic local search methods show a much larger fluctuation of the fitness function than SPORE. Also in this case the quasi-Newton method performs quite poorly, as expected from the simulation results. We should remark that, for each algorithm, the computational time in the real-world experiment was consistent with the values obtained in simulation.

G. Iacca et al. / Applied Soft Computing 23 (2014) 460–473

469

Fig. 6. Simulation results of the second experiment: fitness trend (log scale) of SPORE compared against modern population-based algorithms (top), modern single solution algorithms (center), and classic local search algorithms (bottom).

Finally, Table 3 shows, for each algorithm, the minimum, maximum, mean and median fitness value, as well as its standard deviation over the 300 samples of the considered trajectory. We also indicate the result of the Wilcoxon Rank-Sum test [71] (and its relative p-value), applied with confidence level 0.95 to the pairwise comparison of the fitness trends shown in Fig. 8, taking SPORE

as reference. With reference to the table, the symbols “=” and “+” indicate, respectively, a statistically equivalent performance (i.e., the null hypothesis is accepted) and a better performance of SPORE compared with the algorithm in the row label. It can be seen that the SPORE algorithm proves statistically superior to all the considered population-based algorithms, as well as ISPO, Nelder–Mead,

Table 3 Experimental results of the freehand trajectory depicted in Fig. 3: minimum, maximum, median, and mean fitness value with standard deviation [T2 ] and Wilcoxon Rank-Sum test result (significance level ˛ = 0.05, SPORE taken as reference). Optimizer

Min

Max

Median

Mean

Std. dev.

Wilcoxon

p-value

SPORE MDE-pBX CLPSO CMA-ES ISPO 3SOME nuSA Powell Nelder-Mead Rosenbrock UNCMIN

9.819e−01 4.104e+00 3.697e+01 1.983e+03 6.705e−01 1.624e+00 1.110e+00 1.223e+00 9.489e−01 1.804e+00 5.847e+02

3.015e+06 5.785e+06 6.336e+06 1.457e+07 2.561e+07 2.669e+06 3.090e+06 5.957e+06 8.253e+06 1.008e+07 1.165e+09

4.546e+03 1.283e+04 4.523e+04 6.498e+04 1.842e+04 4.776e+03 5.378e+03 6.593e+03 5.281e+03 1.517e+04 7.703e+06

1.096e+05 6.221e+05 9.081e+05 1.103e+06 7.808e+05 3.076e+05 8.525e+04 2.131e+05 5.498e+05 1.675e+05 2.894e+07

4.374e+05 1.060e+06 1.344e+06 1.942e+06 2.457e+06 7.041e+05 3.896e+05 7.379e+05 1.329e+06 7.603e+05 1.160e+08

+ + + + = = = + + +

2.621e−07 2.985e−17 4.718e−56 3.679e−09 1.822e−01 3.148e−01 1.798e−01 1.024e−02 5.406e−11 2.066e−66

470

G. Iacca et al. / Applied Soft Computing 23 (2014) 460–473

Fig. 7. Simulation results of the first experiment: localization error trend (log scale) of SPORE compared against modern population-based algorithms (top), modern single solution algorithms (center), and classic local search algorithms (bottom).

Rosenbrock and UNCMIN. On the other hand, on this specific trajectory, the SPORE algorithm shows a performance which is statistically equivalent to that of 3SOME, nuSA and Powell. Discussion In both the simulations and the real-world experiments, SPORE proved superior (or at least equivalent), in terms of computational time and accuracy, to all the other algorithms considered for comparison. As mentioned before, the main reason for the greater speed of SPORE has to be found in the number -and type- of computational operations performed by this algorithm, which are much more simple when compared to the operations performed by most modern population-based meta-heuristics or by some classic local-search algorithms (either gradient-based or gradient-free). On the other hand, the better accuracy shown by SPORE confirms its strong exploitation capabilities and indicates that this algorithm might be particularly suitable for solving fitness landscapes

characterized by a single basin of attraction, even in the presence of noise. Also, we expect that the variable-wise perturbations may produce better results on separable problems, although the partial inheritance restart provides a way to perform “diagonal” moves and thus tackle non-separable landscapes. Hence the proposed algorithm might be well-suited in particular for solving efficiently (and quickly) low-dimensional problems, especially with a low level of non-separability and moderate noise. Interestingly, an additional finding of our study is that, for the specific problem of magnetic dipole detection, classic gradient-free optimization methods such as Rosenbrock, Powell and NelderMead algorithms, tend to outperform modern state-of-the-art algorithms, for instance CMA-ES. This result, which might appear counter-intuitive, can be explained in particular considering the relatively low computational budget assigned to the real-time optimization process, in contrast with modern population-based algorithms which usually rely on a large number of fitness evaluations to guarantee high-quality solutions. Finally, the experiment

G. Iacca et al. / Applied Soft Computing 23 (2014) 460–473

471

Fig. 8. Experimental results of the freehand trajectory depicted in Fig. 3: fitness trend (log scale) of SPORE compared against modern population-based algorithms (top), modern single solution algorithms (center), and classic local search algorithms (bottom).

showed that gradient-based quasi-Newton algorithms seem to be not particularly suitable for this kind of problem, most probably because of the noisy features of the fitness landscape. Conclusions and future works In this paper we proposed SPORE, a fast single particle optimization algorithm for the real-time detection of a magnetic moment dipole. Inspired by the search mechanisms of classic single-solution algorithms such as the Hooke–Jeeves pattern search, SPORE performs repeated incrementally smaller variable-wise perturbations in both directions, in the attempt to improve upon the current solution. In addition to that, a simple restart mechanism is used to re-sample the current particle whenever the variable perturbations fail at improving upon it. In order to demonstrate the feasibility of this approach, we implemented a demonstrator setup and we compared the results obtained by SPORE with those of several state-of-the-art and classic alternative methods. Experimental results show that our approach is faster (2.5 ms per sample) and produces on average better results than most of the algorithms under comparison, thus proving a valid

solution for this specific application. The current setup has a worstcase position precision of approximately 6.5 mm and a worst-case orientation precision in the order of tenths of radians. It should be noted, however, that part of this error is due to the noise on the magnetic sensors and to some small errors on the measurement of the sensor positions. In the present study, only the detection and tracking of a single magnetic dipole have been experimentally shown. However, the proposed technique opens up many possibilities for more complex applications. Future studies could continue along two main research lines, which are motion tracking using dipoles and imaging of magnetic objects. As for the first option, a larger number of sensors might be needed to follow multiple dipoles while they move through a confined space. A further development could be using an estimation of the speed of the object to predict the position of the dipole in the future and provide an initial solution in real-time. This system will be most useful in situations where the distance between the dipole and the sensors is not too large, e.g. for pencil tracking. One of the main advantages of a magnetic tracking system would be its ability of detecting objects in the dark or when there are

472

G. Iacca et al. / Applied Soft Computing 23 (2014) 460–473

obstacles. It could even detect dipoles that are inside other (nonmagnetic) objects (e.g. inside the human body, such as a capsule endoscope). The main disadvantage might be instead the relatively short detection range. Regarding the possibility of detecting and imaging arbitrary shaped magnetic objects, the general idea would be to investigate how much information about the shape and size of an object can be extracted based on magnetic field measurements. Possible applications could be the imaging of the magnetic profile of vehicles or other large (magnetic) objects. To obtain a detailed image of such objects, one would need to measure the magnetic field at many different locations. This can be achieved by either a large number of sensors or, alternatively, by scanning the object over the detectable area of an array composed of a few sensors, and recording a series of magnetic field measures while the object moves along the scanning path. Future research could aim to realize an experimental setup to study this magnetic imaging technique. We should note, however, that in both cases a number of assumptions would be needed, since for these complex magnetic objects consisting of multiple magnetic moments with different positions and orientations there could be multiple solutions to the inverse problem. Moreover, due to the higher number of parameters, a larger computation power might also be required. Acknowledgements INCAS3 is co-funded by the Province of Drenthe, the Municipality of Assen, the European Fund for Regional Development and the Ministry of Economic Affairs, Peaks in the Delta. We thank Erik Kallen for the help provided with the hardware setup. References [1] C. Hu, W. Yang, D. Chen, M.Q.-H. Meng, H. Dai, An improved magnetic localization and orientation algorithm for wireless capsule endoscope, in: Annual International Conference of the IEEE Engineering in Medicine and Biology Society, 2008, pp. 2055–2058. [2] C. Hu, M.-H. Meng, M. Mandal, A linear algorithm for tracing magnet position and orientation by using three-axis magnetic sensors, IEEE Trans. Magn. 43 (2007) 4096–4101. [3] K. Levenberg, A method for the solution of certain non-linear problems in least squares, Quart. J. Appl. Math. II (1944) 164–168. [4] W. Yang, C. Hu, M.-H. Meng, S. Song, H. Dai, A six-dimensional magnetic localization algorithm for a rectangular magnet objective based on a particle swarm optimizer, IEEE Trans. Magn. 45 (2009) 3092–3099. [5] D.E. Goldberg, Genetic Algorithms in Search, Optimization and Machine Learning, Addison-Wesley Publishing Co., Reading, MA, USA, 1989. [6] R.C. Eberhart, J. Kennedy, A new optimizer using particle swarm theory, in: Sixth International Symposium on Micromachine and Human Science, 1995, pp. 39–43. [7] C. Cheng, X. Huo, M. Ghovanloo, Towards a magnetic localization system for 3-D tracking of tongue movements in speech-language therapy, in: Annual International Conference of the IEEE Engineering in Medicine and Biology Society, 2009, pp. 563–566. [8] D.R. Jones, C.D. Perttunen, B.E. Stuckman, Lipschitzian optimization without the Lipschitz constant, J. Optim. Theor. Appl. 79 (1993) 157–181. [9] M.J.D. Powell, An efficient method for finding the minimum of a function of several variables without calculating derivatives, Comput. J. 7 (1964) 155–162. [10] A. Nelder, R. Mead, A simplex method for function optimization, Comput. J. 7 (1965) 308–313. [11] C. Xiao, S. Liu, G. Zhou, Real-time localization of a magnetic object with total field data, in: Automation Congress, 2008, pp. 1–4. [12] G. Zhou, C. Xiao, S. Liu, Magnetic moments estimation and localization of a magnetic object based on Hopfield neural network, in: Automation Congress, 2008, pp. 1–5. [13] J.J. Hopfield, Neural networks and physical systems with emergent collective computational abilities, Proc. Natl. Acad. Sci. 79 (1982) 2554–2558. [14] T. Chow, Introduction to Electromagnetic Theory: A Modern Perspective, Jones & Bartlett Learning, MA, USA, 2006. [15] PNI Sensor Corporation, RM3000 & RM2000 Sensor Suite User Manual, 2011, Santa Rosa, CA, USA. [16] F. Camps, S. Harasse, A. Monin, Numerical calibration for 3-axis accelerometers and magnetometers, in: IEEE International Conference on Electro/Information Technology, 2009, pp. 217–221. [17] E. Williams, Geomagnetic Field Calculator, 2012, http://williams.best.vwh.net/ magvar.html

[18] D. Wolff, Using openGL in Java with JOGL, J. Comput. Sci. Colloid 21 (2005) 223–224. [19] J.X. Chen, E.J. Wegman, Foundations of 3D Graphics Programming – Using JOGL and Java3D, Springer-Verlag New York, Inc., Secaucus, NJ, USA, 2006. [20] A. Sokolov, Java-simple-serial-connector, 2011, http://code.google.com/p/ java-simple-serial-connector/ [21] NIST, JAMA: A Java Matrix Package, 1999, http://math.nist.gov/ javanumerics/jama [22] N. Hansen, S.D. Müller, P. Koumoutsakos, Reducing the time complexity of the derandomized evolution strategy with covariance matrix adaptation (CMA-ES), Evol. Comput. 11 (2003) 1–18. [23] C. Igel, T. Suttorp, N. Hansen, A computational efficient covariance matrix update and a (1+1)-CMA for evolution strategies, in: Genetic and Evolutionary Computation Conference, ACM Press, 2006, pp. 453–460. [24] A. Auger, N. Hansen, A restart CMA evolution strategy with increasing population size, in: IEEE Congress on Evolutionary Computation, vol. 2, 2005, pp. 1769–1776. [25] F. Caraffini, G. Iacca, F. Neri, L. Picinali, E. Mininno, A CMA-ES super-fit scheme for the re-sampled inheritance search, in: IEEE Congress on Evolutionary Computation, 2013, pp. 1123–1130. [26] I. Loshchilov, CMA-ES with restarts for solving CEC 2013 benchmark problems, in: IEEE Congress on Evolutionary Computation, 2013, pp. 369–376. [27] R. Mallipeddi, P.N. Suganthan, Q.K. Pan, M.F. Tasgetiren, Differential evolution algorithm with ensemble of parameters and mutation strategies, Appl. Soft Comput. 11 (2011) 1679–1696, The Impact of Soft Computing for the Progress of Artificial Intelligence. ´ M. Mernik, V. Zˇ umer, Self-adapting control [28] J. Brest, S. Greiner, B. Boˇskovic, parameters in differential evolution: a comparative study on numerical benchmark problems, IEEE Trans. Evol. Comput. 10 (2006) 646–657. [29] R. Mallipeddi, G. Iacca, P.N. Suganthan, F. Neri, E. Mininno, Ensemble strategies in compact differential evolution, in: IEEE Congress on Evolutionary Computation, 2011, pp. 1972–1977. [30] M.G.H. Omran, A. Salman, A.P. Engelbrecht, Self-adaptive differential evolution, in: Computational Intelligence and Security, volume 3801 of Lecture Notes in Computer Science, Springer, 2005, pp. 192–199. [31] K.E. Parsopoulos, Cooperative micro-differential evolution for highdimensional problems, in: Genetic and Evolutionary Computation Conference, 2009, pp. 531–538. [32] A.K. Qin, V.L. Huang, P.N. Suganthan, Differential evolution algorithm with strategy adaptation for global numerical optimization, IEEE Trans. Evol. Comput. 13 (2009) 398–417. [33] A.K. Qin, P.N. Suganthan, Self-adaptive differential evolution algorithm for numerical optimization, in: IEEE Congress on Evolutionary Computation, volume 2, 2005, pp. 1785–1791. [34] F. Caraffini, F. Neri, J. Cheng, G. Zhang, L. Picinali, G. Iacca, E. Mininno, Super-fit multicriteria adaptive differential evolution, in: IEEE Congress on Evolutionary Computation, 2013, pp. 1678–1685. [35] S. Rahnamayan, H.R. Tizhoosh, M.M. Salama, Opposition-based differential evolution, IEEE Trans. Evol. Comput. 12 (2008) 64–79. [36] X. Li, X. Yao, Cooperatively coevolving particle swarms for large scale optimization, IEEE Trans. Evol. Comput. 16 (2012) 210–224. [37] J.J. Liang, A.K. Qin, P.N. Suganthan, S. Baskar, Comprehensive learning particle swarm optimizer for global optimization of multimodal functions, IEEE Trans. Evol. Comput. 10 (2006) 281–295. [38] Q. Yuan, F. Qian, W. Du, A hybrid genetic algorithm with the Baldwin effect, Inform. Sci. 180 (2010) 640–652. [39] S. Ghosh, S. Das, D. Kundu, K. Suresh, A. Abraham, Inter-particle communication and search-dynamics of lbest particle swarm optimizers: an analysis, Inform. Sci. 182 (2012) 156–168. [40] M.A. Montes de Oca, T. Stutzle, M. Birattari, M. Dorigo, Frankenstein’s PSO: a composite particle swarm optimization algorithm, IEEE Trans. Evol. Comput. 13 (2009) 1120–1132. [41] Y. Shi, H. Liu, L. Gao, G. Zhang, Cellular particle swarm optimization, Inform. Sci. 181 (2011) 4460–4493 (Special Issue on Interpretable Fuzzy Systems). [42] H. Wang, Z. Wu, S. Rahnamayan, Y. Liu, M. Ventresca, Enhancing particle swarm optimization using generalized opposition-based learning, Inform. Sci. 181 (2011) 4699–4714. [43] Y. Wang, B. Li, T. Weise, J. Wang, B. Yuan, Q. Tian, Self-adaptive learning based particle swarm optimization, Inform. Sci. 181 (2011) 4515–4538. [44] W.-N. Chen, J. Zhang, Y. Lin, N. Chen, Z.-H. Zhan, H.-H. Chung, Y. Li, Y.-H. Shi, Particle swarm optimization with an aging leader and challengers, IEEE Trans. Evol. Comput. 17 (2013) 241–258. [45] N. Noman, H. Iba, Accelerating differential evolution using an adaptive local search, IEEE Trans. Evol. Comput. 12 (2008) 107–125. [46] D.H. Kim, A. Abraham, J.H. Cho, A hybrid genetic algorithm and bacterial foraging approach for global optimization, Inform. Sci. 177 (2007) 3918– 3937. [47] C.W. Ahn, J. An, J.-C. Yoo, Estimation of particle swarm distribution algorithms: Combining the benefits of PSO and EDAs, Inform. Sci. 192 (2012) 109– 119. [48] M.N. Le, Y.S. Ong, Y. Jin, B. Sendhoff, Lamarckian memetic algorithms: local optimum and connectivity structure analysis, Memetic Comput. J. 1 (2009) 175–190. [49] R. Meuth, M.H. Lim, Y.S. Ong, D.C. Wunsch-II, A proposition on memes and meta-memes in computing for higher-order learning, Memetic Comput. J. 1 (2009) 85–100.

G. Iacca et al. / Applied Soft Computing 23 (2014) 460–473 [50] D. Molina, M. Lozano, C. Garcia-Martinez, F. Herrera, Memetic algorithms for continuous optimization based on local search chains, Evol. Comput. 18 (2010) 27–63. [51] G. Iacca, F. Caraffini, F. Neri, E. Mininno, Single particle algorithms for continuous optimization, in: IEEE Congress on Evolutionary Computation, 2013, pp. 1610–1617. [52] Q.C. Nguyen, Y.S. Ong, M.H. Lim, A probabilistic memetic framework, IEEE Trans. Evol. Comput. 13 (2009) 604–623. [53] F. Neri, E. Mininno, Memetic compact differential evolution for cartesian robot control, IEEE Comput. Intell. Mag. 5 (2010) 54–65. [54] F. Neri, G. Iacca, E. Mininno, Disturbed exploitation compact differential evolution for limited memory optimization problems, Inform. Sci. 181 (2011) 2469–2487. [55] G. Iacca, F. Caraffini, F. Neri, Memory-saving memetic computing for pathfollowing mobile robots, Appl. Soft Comput. 13 (2013) 2003–2016. [56] F. Caraffini, F. Neri, G. Iacca, A. Mol, Parallel memetic structures, Inform. Sci. 227 (2013) 60–82. [57] F. Neri, E. Mininno, G. Iacca, Compact particle swarm optimization, Inform. Sci. 239 (2013) 96–121. [58] G. Iacca, Distributed optimization in wireless sensor networks: an island-model framework, Soft Comput. (2013) 1–21. [59] G. Iacca, F. Caraffini, F. Neri, Multi-strategy coevolving aging particle optimization, Int. J. Neural Syst. 24 (2014) 1–19. [60] F. Caraffini, F. Neri, B. Passow, G. Iacca, Re-sampled inheritance search: high performance despite the simplicity, Soft Comput. (2013) 1–22.

473

[61] K.V. Price, R. Storn, J. Lampinen, Differential Evolution: A Practical Approach to Global Optimization, Springer-Verlag New York, Inc., Secaucus, NJ, USA, 2005. [62] R. Hooke, T.A. Jeeves, Direct search solution of numerical and statistical problems, J. ACM 8 (1961) 212–229. [63] S. Islam, S. Das, S. Ghosh, S. Roy, P. Suganthan, An adaptive differential evolution algorithm with novel mutation and crossover strategies for global numerical optimization, IEEE Trans. Syst. Man Cybern. B: Cybern. 42 (2012) 482–500. [64] N. Hansen, The CMA Evolution Strategy, 2012, http://www.lri.fr/hansen/ cmaesintro.html [65] G. Iacca, F. Neri, E. Mininno, Y.S. Ong, M.H. Lim, Ockham’s razor in memetic computing: three stage optimal memetic exploration, Inform. Sci. 188 (2012) 17–43. [66] J. Zhou, Z. Ji, L. Shen, Simplified intelligence single particle optimization based neural network for digit recognition, in: Chinese Conference on Pattern Recognition, 2008, pp. 1–5. [67] Z. Xinchao, Simulated annealing algorithm with adaptive neighborhood, Appl. Soft Comput. 11 (2011) 1827–1836. [68] R.P. Brent, Algorithms for minimization without derivatives, Prentice-Hall, Englewood Cliffs, NJ, 1973. [69] H.H. Rosenbrock, An automatic method for finding the greatest or least value of a function, Comput. J. 3 (1960) 175–184. [70] R.B. Schnabel, J.E. Koonatz, B.E. Weiss, A modular system of algorithms for unconstrained minimization, ACM Trans. Math. Softw. 11 (1985) 419–440. [71] F. Wilcoxon, Individual comparisons by ranking methods, Biomet. Bull. 1 (1945) 80–83.