An agent-Oriented Hierarchic Strategy for Solving Inverse Problems

Int. J. Appl. Math. Comput. Sci., 2015, Vol. 25, No. 3, 483–498 DOI: 10.1515/amcs-2015-0036 AN AGENT–ORIENTED HIERARCHIC STRATEGY FOR SOLVING INVERSE PROBLEMS MACIEJ SMOŁKAa,∗,ROBERT SCHAEFERa,MACIEJ PASZY ´ NSKIa,DAVID PARDOb,c,d, JULEN ´ ALVAREZ-ARAMBERRIb,e aDepartment of Computer Science AGH University of Science and Technology, Al. Mickiewicza 30, 30-059 Krak´ow, Poland e-mail: {smolka,schaefer,paszynsk}@agh.edu.pl bDepartment of Applied Mathematics, Statistics, and Operational Research University of the Basque Country (UPV/EHU), Barrio Sarriena S/N, 48940 Leioa (Bizkaia), Spain e-mail: {dzubiaur,julen.alvarez.aramberri}@gmail.com cBasque Center for Applied Mathematics (BCAM), Alameda de Mazarredo 14, 48009 Bilbao, Spain dIkerbasque (Basque Foundation for Sciences), Maria Diaz de Haro 3, 48013 Bilbao, Spain eLaboratory of Mathematics and Their Applications University of Pau (UPPA), Avenue de l’Universit´e BP 576, 64012 Pau cedex, France The paper discusses the complex, agent-oriented hierarchic memetic strategy (HMS) dedicated to solving inverse parametric problems. The strategy goes beyond the idea of two-phase global optimization algorithms. The global search performed by a tree of dependent demes is dynamically alternated with local, steepest descent searches. The strategy offers exceptionally low computational costs, mainly because the direct solver accuracy (performed by the hp-adaptive finite element method) is dynamically adjusted for each inverse search step. The computational cost is further decreased by the strategy employed for solution inter-processing and fitness deterioration. The HMS efficiency is compared with the results of a standard evolutionary technique, as well as with the multi-start strategy on benchmarks that exhibit typical inverse problems’ difficulties. Finally, an HMS application to a real-life engineering problem leading to the identification of oil deposits by inverting magnetotelluric measurements is presented. The HMS applicability to the inversion of magnetotelluric data is also mathematically verified. Keywords: inverse problems, hybrid optimization methods, memetic algorithms, multi-agent systems, magnetotelluric data inversion. 1. Introduction Inverse problems form an important area of contemporary research related to fundamental problems in science and engineering. Among its numerous applications one can find non-invasive geophysical exploration, tumor characterization, and analysis of unknown materials. Parametric inverse problems are usually formulated as optimization ones, where the objective is to minimize the misfit between the simulated and measured forward solutions. When solving such problems, one usually ∗Corresponding author faces some significant obstacles such as ill-conditioning, existence of multiple local minima (multi-modality), and possibly low regularity of the misfit functional. All of them significantly reduce the usefulness of convex optimization methods (such as gradient-based ones), as well as simple stochastic mechanisms (Monte Carlo and simple evolutionary schemes) because of •the lack of guarantee of finding all solutions, •enormous computational cost. The main goal of this work is to obtain a global inverse solver capable of finding many global and/or local 484 M. Smołka et al. solutions (misfit minima) with a satisfactory accuracy and an acceptable computational cost. In order to achieve this goal, the paper combines several different ideas: the hierarchic genetic strategy (HGS) (see, e.g., Schaefer and Kołodziej, 2003; Wierzba et al., 2003) to decrease the cost of the main part of a global search (which consists of finding the basins of attraction), cluster-based fitness deterioration (see, e.g., Beasley et al., 1993; Obuchowicz, 1997; Wolny and Schaefer, 2011), memetic algorithms (see, e.g., Neri et al., 2012) composing various techniques into a single population-based stochastic strategy in order to gain efficiency and flexibility, and the evolutionary multi-agent system (EMAS) (see, e.g., Cetnarowicz et al., 1996), allowing more flexible, asynchronous processing and easy embedding of local convex searches. The work presented in this paper is an extension of the hybrid model involving the global genetic search phase followed by the local gradient phase, as applied for the inverse resistivity logging measurement simulations with direct current electrodes (Gajda-Zag´orska et al., 2015), as well as for the identification of the Young modulus for imprint nanolithography (Barabasz et al., 2014). The proposed strategy develops dynamically a tree of dependent populations (demes) searching with various levels of accuracy that grow from the root to the leaves. Two types of individuals are utilized: passive individuals, containing candidate solutions only, and active individuals, consisting of computational software agents. Both active and passive individuals are gathered in demes governed by structural agents assigned to the nodes of the population tree. Demes of passive individuals are traditionally evolved by using common selection and genetic operators. Active individuals (agents) compete inside a deme one with another to perform their actions producing offspring agents with genotypes determined by traditional genetic operators or by a convex optimization process. The search accuracy is that associated to the solution of direct problems using an hp-adaptive finite element method (hp-FEM), where h (height) refers to the element size, and pto the polynomial order of approximation, which can be adapted/modified throughout the computational grid (Demkowicz, 2006). The low computational cost is obtained primarily by an economical global search performed by the tree of demes, because the more detailed local searches are activated mainly in the promising regions found by the parental populations. The total number of individuals is significantly lower than in the case of the traditional, single population evolutionary search. The main computational cost decrement is achieved via the common inverse and forward error scaling, i.e., the rough global search is performed with a low misfit accuracy whereas the accuracy increases in the refined searches performed by branches and leaves. Such a policy minimizes the number of expensive hp-FEM solver calls, which are necessary for accurate misfit evaluation. Finally, the agent-oriented architecture based on the EMAS idea allows economic and flexible invocation of local gradient methods at least once for a single solution and only if they outperform the stochastic search. Moreover, the agent-based architecture facilitates parallel evolution of deme populations. We have verified our strategy by a series of benchmark problems that mimic the optimization landscape obtained in real-life inverse problems. The first phase of benchmark tests leads to establishing the set of the strategy’s parameters that assure its most economic operation. In the second phase, we compare results of our proposed strategy with those of other global optimization methods executed using a comparable budget (computer resources). In this paper, we compare the HMS with two standard global optimization methods, i.e., the simple evolutionary algorithm (SEA) and multi-start. We have already performed a benchmark-based comparison with a version of the HGS (Smołka and Schaefer, 2014). We conclude the paper with the application of our agent-based system to solve an inverse problem consisting in identifying the resistivity of the underground formation by using magnetotelluric measurements. The magnetotelluric (MT) method is a passive electromagnetic (EM) exploration technique which allows us to determine the resistivity distribution in the subsurface of the area of interest on scales varying from a few meters to hundreds of kilometers. Commercial applications include hydrocarbon, geothermal, underground water monitoring, and, more recently, monitoring CO2sequestration in the subsurface. The method does not require artificial power sources. Instead, natural power sources are used to induce the phenomena. These natural sources are nothing but electric currents within the ionosphere created by deformations of the magnetosphere. The motivation for the choice of the MT problem was the detection of several local minima in an MT inversion performed by means of a classical gradient method ( ´ Alvarez-Aramberri et al., 2013). The results of all comparative tests show a much higher quality and number of extremes found by the proposed agent-based system than by other reference methods. Single population evolutionary algorithms are able to find only a single solution within the assumed budget. A multi-start method, which consists in parallel execution of local gradient-type methods from a number of random starting points, cannot in general find satisfactory solutions, because within the assumed budget only a small number of local processes can be invoked. The paper concentrates on the algorithmic aspects of the HMS, so we do not discuss implementation-related issues here. Some notes on a sample implementation can be found in the work of Smołka and Schaefer (2014) and we refer the interested reader there. An agent-oriented hierarchic strategy for solving inverse problems 485 2. HMS architecture The main idea of the HMS is to provide a global optimization tool especially suited for solving difficult inverse problems. The considered problems are difficult because of their inherent multi-modality accompanied by the nontrivial computational cost of direct problem solution, which is necessary for evaluating the objective that is the misfit between the observed and computed values of a quantity of interest. Nevertheless, inverse problems have also some advantageous features. First, their theoretical global minimum value is well known (and equal to zero). Although in practice this theoretical value is never attained because of noise, modeling errors, etc., we know at least that the objective function is bounded from below by zero and that its takes values close to this bound. This knowledge can be used, e.g., in the construction of stopping conditions for stochastic evolution. Second, in some importantcases, the cost of the direct problem solution can be modulated by an assumed accuracy of the solution: it is the case of hp-FEM direct solvers (Demkowicz et al., 2007). As a global optimization tool, the HMS tries to combine the high-level exploratory ability with the accuracy and efficiency of a local optimization method. In contrast to two-phase methods in which local searches are executed right after the completion of the global phase, the HMS follows the overall idea of memetic algorithms, i.e., it intermixes local-optimization-oriented mechanisms into a global stochastic search machinery. The global part follows the multi-population evolutionary approach introduced by the HGS (Schaefer and Kołodziej, 2003). Namely, the global search is performed by a collection of genetic populations. The populations can evolve in parallel, but they are not mutually independent. The structure of the dependency relation is hierarchical (i.e., tree-like; see Fig. 1) with a restricted number of levels. The HGS proved to have Level 1 Level 2 Level 3 root deme branch demes leaf demes U1 genetic spaces low accuracy high accuracy U2 U3 Fig. 1. HGS-like evolutionary population tree. considerable exploratory capabilities together with a good search accuracy, especially with floating-point phenotype encoding (Wierzba et al., 2003). The HMS naturally tries to retain these abilities while going beyond the HGS in some aspects. First of all, it adds local optimization to the set of operations applied to the genetic individuals. But this is performed carefully in order to avoid the premature population convergence on the one hand and high cost of running instances of a local method from inappropriate points on the other. Namely, some genetic individuals (but not necessarily all of them) receive an identity and some intelligence, hence becoming independent agents in a multi-agent system (MAS). The decision of performing the local search becomes their own responsibility. Moreover, the demes are managed by special controller agents. The idea of turning a passive genetic individual into an intelligent agent has some further consequences. We have to redefine the genetic operations in such a way that they can be applied to agents. While it is straightforward for the mutation and the crossover (although in this case a new agent is activated), the selection cannot be performed in the simple genetic (or evolutionary) way. Instead, we follow the lines of the EMAS (Cetnarowicz et al., 1996; Byrski et al., 2013), thus performing an operation analogous to the tournament selection but realized as a two-agent rendezvous. In the sequel, we shall present the structure of the HMS starting with a description of HMS agent types. 2.1. HMS agent types. The main HMS agent types along with their interrelations are shown in Fig. 2. The Master Agent (MA) is a global system coordinator. Deme Agents (DAs) manage evolutionary populations at various levels of the deme tree. DA specializations differ primarily in the type of population they can own. Evolutionary Agents (EAs) hold simple passive collections of chromosomes, whereas Local Agents (LAs) coordinate groups of Computational Agents (CAs). Any CA, apart from holding an immutable genotype, can perform one of the available actions, such as local optimization method execution. The primary responsibility of the Objective Agent (OA) is objective computation, which typically involves calling an external direct problem solver. In the following, we provide detailed descriptions of HMS agent types. MasterAgent ObjectiveAgent cache DemeAgent accuracyLevel EvolutionaryAgent passivePopulation LocalAgent ComputationalAgent genotype {readOnly} lifeEnergy * 1 * ExternalSolver Fig. 2. HMS agent types (UML class diagram). 486 M. Smołka et al. Master Agent (MA). As a global system coordinator, it is started as a first agent in the HMS MAS. Its responsibilities include performing system initialization such as the activation of other essential agents, i.e., the Objective Agent and a Deme Agent of the deme-tree root. After the initialization, the Master Agent starts the global loop of deme coordination and checks if the global stopping condition is satisfied. The deme coordination follows the lines of FIPA Contract-Net (FIPA, 2002). It begins by sending a call for proposals (CFP) to all active Deme Agents. Then, the MA waits for DA proposals and accepts those that are not in conflict. This is shown in Algorithm 1. DA proposals are in conflict when their Algorithm 1. Master Agent (MA) algorithm. 1: create OA 2: create root location DA 3: repeat 4: send CFP to all active DAs 5: receive proposals from DAs 6: accept all non conflicting proposals 7: request fitness deterioration from OA 8: until global stop condition is satisfied. corresponding activities cannot be executed in parallel. Using the terminology of Byrski et al. (2013), they are not local. In the current HMS realization, all actions are local (note that, in contrast to Byrski et al. (2013), we do not consider the migration), so all DA proposals can be accepted by the MA. Deme Agent (DA). It is a deme-tree node coordinator. Each deme has an associated level of computational accuracy stored as a property of the corresponding Deme Agent. In fact, the Deme Agent is an abstract class with two different specializations: the Evolutionary Agent and the Local Agent. Evolutionary Agent (EA). This is a simple (passive) evolutionary population owner. Periodically, after receiving the permission from the Master Agent, it evolves its population for a fixed number of generations (this sequence of genetic epochs is called a metaepoch), and then sprouts a new deme from the current best individual unless the sprout condition is not satisfied (see Algorithm 2). Note that similar agents form the structure of the globally balanced HGS (Jojczyk and Schaefer, 2009). Creating the initial population in Line 2 has two different meanings. If an EA contains the tree root deme, then the initial population is sampled using the uniform probability distribution. Otherwise, the EA is itself sprouted using its parent’s best individual as a seed. In such a case, the initial population is sampled using the Algorithm 2. Evolutionary Agent (EA) algorithm. 1: set accuracy level 2: create initial deme population 3: repeat 4: send proposal to MA 5: if MA has accepted proposal then 6: for all epochs in metaepoch do 7: perform selection 8: perform crossover and mutation 9: for all created individual do 10: request objective computation from OA with stored accuracy 11: end for 12: end for 13: if best individual satisfies sprout condition then 14: sprout new child DA from best individual 15: end if 16: end if 17: until local stop condition is satisfied normal distribution centered in the seed individual with the standard deviation depending upon the tree level. Local Agent (LA). The Local Agent owns a population of Computational Agents and acts as a local scheduler of their actions. Namely, it receives action proposals from Computational Agents, selects one of them according to a probability distribution, sends a proposal to the Master Agent and, if the proposal is accepted, lets the selected Computational Agent perform its action (see Algorithm 3). The Local Agent’s responsibilities include also some action coordination, such as checking if a selected sprout action is allowed. Creating the initial population is performed analogously as in the case of an Evolutionary Agent, i.e., by sampling using a proper probability distribution (different for root and non-root demes). The only difference is the type created individuals: passive chromosomes in the case of an Evolutionary Agent and active Computational Agents in the case of Local Agents. Computational Agent (CA). It is an active individual of the HMS genetic population. It owns an immutable genotype consisting of an encoded domain point (a chromosome) and a level of computational precision. The precision level must be consistent with the owning Local Agent’s level. The mutable part of a Computational Agent’s state includes a nonnegative memetic parameter called life energy. It is exchanged during the ComputationalAgent’s action execution such that the total energy remains constant within each deme. Only agents with positive life energy are considered active (alive) and take part in system evolution. There exists a set of actions An agent-oriented hierarchic strategy for solving inverse problems 487 Algorithm 3. Local Agent (LA) algorithm. 1: set accuracy level 2: create initial deme population 3: repeat 4: send CFP to all active CAs 5: receive action proposals from CAs and choose one 6: send corresponding proposal to MA 7: if MA has accepted the proposal then 8: if CA action creates new individual then 9: create new CA 10: else if chosen action is SPROUT then 11: if sprouting can be performed then 12: create new child DA 13: end if 14: end if 15: end if 16: until local stop condition is satisfied from which an active Computational Agent chooses one at a time to perform. Which actions are actually available depends on parameters such as life energy, the objective value, the precision level, and others. Finally, the action is performed only if permitted by the owning Local Agent (see Algorithm 4). Algorithm 4. Computational Agent (CA) algorithm. 1: set accuracy level according to LA’s level 2: request objective computation from OA with stored accuracy 3: repeat 4: remain temporarily inactive 5: until objective value is obtained 6: while life energy >0do 7: receive CFP from owning LA 8: choose an available action 9: send the corresponding proposal to LA 10: if received permission from LA then 11: perform chosen action 12: update life energy 13: end if 14: end while{CA permanently inactive} There is an energy quantum related to each action, which is spent (during GET it can sometimes be gained) by a Computational Agent during action execution. Currently, the following actions are considered (cf. Byrski et al., 2013): GET, MUTATE, CROSSOVER, LOCOPT and SPROUT. In all tests prepared for this paper, we set all action energies to 1and the initial CA energy to 5. The GET action is a two-agent stochastic duel during which proper quantum energy moves from the loser to the winner. A Computational Agent with a lower (i.e., closer to the global minimum) objective value has more chances to win. MUTATE and CROSSOVER are straightforward counterparts of the corresponding genetic (or evolutionary) operations such as, e.g., the normal mutation and the arithmetic crossover. The SPROUT action is inspired by the child branch sprouting operation, which is fundamentalin the HGS (Schaefer and Kołodziej, 2003). In the HMS, it produces a new deme together with its Deme Agent and an initial population of Computational Agents. Obviously, SPROUT makes no sense at the leaf level, where it can be optionally replaced with LOCOPT. LOCOPT is local optimization method execution started from the agent’s decoded chromosome. In the current realization, LOCOPT is allowed only at the leaves. Action selection is determined by a given probability distribution. The probability of LOCOPT can be computed using a formula like the following: pLOCOPT =p0+(1−p0)s 1+objective,(1) where 0≤p0<1is the guaranteed probability and s>0is close but not equal to 1to allow selecting other actions. Thus, the better the objective value, the more strongly LOCOPT is preferred. Other actions available according to the current life energy receive equal remaining probability. The same formula can also be used for computing SPROUT probability at non-leaf levels. Objective Agent (OA). In a real HMS application (i.e., in solving inverse problems), the objective value is computed externally by a specialized direct solver. The responsibility of an Objective Agent (typically, one in the whole system) is to provide a proper solver gateway, i.e., to execute the solver process (or several parallel processes) properly and to transfer the input data to the solver and the output back to the HMS. Of special interest for us is the case of computing the objective by means of a direct hp-FEM solver when the direct problem solution is Lipschitz continuous with respect to the parameters. This property is not straightforward and has to be proved in any particular case. In Section 4.4 this is shown for the magnetotelluric problem (see Remark 1), which is our real-world test case. Then, we can adapt the solver accuracy to the assumed accuracy of HMS tree demes; see Algorithm 5. Finally, we know the dependency between the solver accuracy and the computational cost of the direct problem solution (for details, see Barabasz et al., 2014; Gajda-Zag´orska et al., 2015) which in turn is the main unit component of the overall HMS cost. Therefore, we can optimize the overallcost by modulating the deme accuracy. Additional Objective Agent activities may include caching solver results, solver instance pooling (in the case of parallel execution) and scheduling objective computations according to a sophisticated optimizing 488 M. Smołka et al. Algorithm 5. Objective computation with hp-FEM direct solver (see Section 4.4). 1: input: tree level j, parameter value 2: compute relative FEM error erel 3: while erel <Ratio(j)do 4: execute 1step of hp adaptation 5: solve problem on new fine and coarse meshes 6: computenewvalueoferel 7: end while 8: return objective computed by means of final mesh policy (e.g., a diffusion-based one (Grochowski et al., 2006)). Quite a special kind of additional OA activity is the deterioration of the objective function. The latter is a proper objective modification that leads to the leveling of central regions of attraction basins of already found local and global minima (Beasley et al., 1993; Obuchowicz, 1997). The idea is to discourage the evolutionary individuals from exploring already well-recognized areas. In this paper, we adopt the cluster-based fitness deterioration technique described by Wolny and Schaefer (2011). Namely, after a request from the MA, all gathered objective data (in this case, the points at which the objective has been computed) at selected accuracy levels are clustered (currently, using the DBSCAN algorithm (Ester et al., 1996)). Afterwards, for each recognized cluster, we construct the ellipsoidal cluster extension CE ={x∈RN:(x−x)TΣ−1(x−x)≤1} determined by its center xand a symmetric matrix Σ, which in our case is the cluster unbiased sample covariance matrix (cf. Wolny and Schaefer, 2011). Then, each cluster extension contributes to the deterioration by the following formula: f(x):=f(x)+fmaxΨ(x),(2) where Ψis a well-known bump function Ψ(x) =⎧ ⎨ ⎩ exp 1−1 1−(x−x)TΣ−1(x−x)for x∈CE, 0otherwise. (3) and fmax is the maximal objective value over the cluster. Ψis an infinitely smooth function and its support is exactly the considered ellipsoid. Note that, in contrast to the standard (cf. Wolny and Schaefer, 2011), we do not make use of simple Gaussian functions because they are positive everywhere, which destroys the desired zero-value feature of inverse problem global minima. 2.2. Population structure. As stated before, the HMS genetic population is decomposed into dependent demes forming a dynamically changing tree of the fixed maximal depth m. Genetic individuals, i.e., computing agents, located at the tree levels close to the root, perform the chaotic and inaccurate search. When approaching the leaves, the search becomes more and more focused and the accuracy is increased (see Fig. 1). The variability of the search accuracy results from the diversity of the genotype encoding precision used at different tree levels. The latter depends on the encoding type. In the case of binary encoding (as in the simple genetic algorithm), it can be achieved by the binary genotype length variation, whereas in the case of real number encoding (as in the simple evolutionary algorithm), it can be realized by appropriate phenotype scaling. The latter case is used in the prototype implementation of the HMS, so here we present some details. The description follows those presented in existing papers (Wierzba et al., 2003; Jojczyk and Schaefer, 2009). In real number encoding, both phenotypes and genotypes are vectors from RN. We assume that the solution domain is a box D=[a1,b 1]×···×[aN,b N], and we take a sequence of scaling factors ηi∈Rsuch that η1>η 2> ...η m−1>η m=1. Then, the genetic universum at the tree level jis Uj=0,b1−a1 ηj×···×0,bN−aN ηj,(4) and the encoding mapping at the level jis defined as Dx−→ xk−ak ηjN k=1 ∈Uj.(5) Moreover, we define the scaling mapping, scalei,j :Uix→ ηi ηj x∈Uj. In such genetic universa, the search at lower levels is more chaotic because the mutation acts stronger, and less precise because of limitations in the real number representation. It is possible to use various genetic operators in such encoding. Among the most important ones, we select the normal mutation yi=xi+N(0,σmut j),i=1,...,N, where N(0,σmut j)is a normally distributed random variable with standard deviation σmut jset separately for each level j, and the arithmetic crossover yi=x1 i+U([0,1])(x2 i−x1 i),i=1,...,N, An agent-oriented hierarchic strategy for solving inverse problems 489 where U([0,1]) is a random variable distributed uniformly over the interval [0,1]. Both operators are used in our sample implementation. Furthermore, we exploit the classical fitness-proportional (roulette-wheel) selection in passive populations(on Evolutionary Agents), additionally preserving the best individual of each generation. A newly sprouted deme’s population is sampled according to the N-dimensional Gaussian distribution centered at the properly encoded fittest individual of the parent process with the diagonal covariance matrix taking values (σsprout j)2on the diagonal. The sprout cannot be performed in population Pat level j if there exists a population Pat level j+1such that |y−scalei,i+1(y)|<c j,(6) where yis the best individual in P,yis the average phenotype of P,andcjis a branch comparison constant. Finally, it should be mentioned that further utilization of the knowledge gathered during multi-level enhanced genetic evolution is possible by means of the clustering technique, in which better approximation of attraction basins of the local minima can be developed, allowing yet more precise application of local optimization methods. 3. Benchmark tests In order to prove HMS abilities to find the global minimum in multi-modal cases, we performed two benchmark tests. The test setup in both the cases was as follows. The HMS tree had three levels with Evolutionary Agents at the root and middle levels, and Local Agents at the leaf level, which seems to be quite a standard layout for an overall optimization. The proportional (roulette-wheel) selection, the normal mutation and the arithmetic crossover were chosen as genetic operators. First, we ran the HMS against both functions 30 times till the global minimum was reached with the assumed accuracy of 10−4. The choice of such a type of test influenced the setting of the HMS stopping conditions. Namely, the global stopping condition was satisfied if a leaf approached the global minimum with the given accuracy, whereas a leaf stopping condition was satisfied if the leaf approached the global minimum or if a fixed number of its consecutive metaepochs were ineffective, i.e., there was no significant variation in the leaf’s population average fitness. To make this stopping condition applicable to active populations, we need to adapt the notion of the genetic epoch. Namely, in the case of an LA population, it is simply a CA action execution sequence of the length equal to the initial population size. After the HMS terminated, we counted its objective calls, which set the computational budget for two comparative classical stochastic optimization methods: the simple evolutionary algorithm (SEA) and the multi-start. The comparative methods were also run 30 times. The first benchmark was the 10-dimensional Ackley path function with domain [−5,5]10. It is a standard global optimization test with one hard-to-find global zero-valued minimum surrounded by numerous other local minima with greater values. This makes it similar to a few important inverse problem objectives, such as the magnetotelluric problem (see Section 4). A special feature of the Ackley function is its flatness outside the narrow attraction basin of the global minimum, which makes it still more similar to the MT problem. The execution parameters for the 10D Ackley function are summarized in Table 1. Note that the metaepoch length parameter is specially adapted to Local Agents (see above). Similarly, the population size in this case is not constant; in our simulations it varied between 4and 8. Table 1. HMS execution parameters (Ackley 10D). Root Middle Leaf Population (initial) 50 10 4 Metaepoch length 2 2 2 Encoding scale 4.0 2.0 1.0 Mutation rate 0.2 0.05 0.01 Crossover rate 0.5 0.5 0.5 Mutation std. dev. 5.0 1.0 0.2 Sprout std. dev. –2.0 1.0 Sprout min. dist. –2.0 1.0 The obtained HMS objective call means are shown in Table 2. The cost of local method application is included in the leaf level cost. The averages from Table 2 Table 2. Average number of objective evaluations (Ackley 10D). Root Middle Leaf Total 1022.8584.74493.66101.1 allowed us to compute the predicted cost for the SEA and multi-start. Taking the average cost of running the local method from HMS executions, we set the starting pool size for multi-start to 70. Similarly, estimating the average cost of running the evolutionary algorithm, also on the basis of HMS runs, we set the initial SEA population size to 100. The SEA genetic parameters were set according to the corresponding HMS parameters from the leaf level, i.e., the crossover rate was set to 0.5, the mutation rate to 0.01, and the mutation standard deviation to 0.2.Table3 shows the average obtained minimum values for all the three methods. Neither the SEA nor the multi-start ever succeeded in reaching the actual global minimum, which the HMS managed to do in every run. Full statistics are shown in a concise way in Fig. 3. Horizontal bars indicate minimum, mean and maximum, respectively. The width of the ‘violin’ indicates the distribution of values. The second benchmark was the product of three reflected and vertically translated Gaussian functions over 490 M. Smołka et al. Table 3. Average values of computed global minimum (Ackley 10D). HMS Multi-start SEA Average 5.27 ·10−93.17 2.57 Best 1.96 ·10−10 1.65 1.3 Fig. 3. Global minima statistics (Ackley 10D). the 4-dimensional box [−5,5]4, f2(x)= 3  i=1 1−exp −(x−mi)TAi(x−mi),(7) where m1=(3,2,1,0),m2=(1,−3,−1,3),m3= (−2.5,2,2,2) and A1=⎡ ⎢ ⎢ ⎣ 1000 0100 0010 0001 ⎤ ⎥ ⎥ ⎦ , A2=⎡ ⎢ ⎢ ⎣ 0.2000 00.10 0 000.10 0000.2 ⎤ ⎥ ⎥ ⎦ , A3=⎡ ⎢ ⎢ ⎣ 1.50 0 0 0200 001.50 0001 ⎤ ⎥ ⎥ ⎦ . Benchmark f2has three separate zero-valued global minima m1,m2and m3with no other local minima. The difficulty of finding the global minima is graded. Here m2has the broadest attraction basin, so it is fairly easily reachable. Then m1is steeper, hence more difficult to find. The main problem is reaching m3. Outside their attraction basins, one can find quite large plateaus, which can cause serious trouble for gradient optimization methods. Thus, the aim of this test is to check the compared methods’ abilities of finding all global minima. The HMS execution parameters for benchmark f2 are gathered in Table 4. The obtained HMS objective call Table 4. HMS execution parameters (three Gaussians). Root Middle Leaf Population (initial) 50 10 4 Metaepoch length 2 2 2 Encoding scale 4.0 2.0 1.0 Mutation rate 0.4 0.1 0.01 Crossover rate 0.5 0.5 0.5 Mutation std. dev. 10.0 1.0 0.1 Sprout std. dev. –2.0 1.0 Sprout min. dist. –2.0 1.0 means are shown in Table 5. This time, to make things Table 5. Average number of objective evaluations (three Gaussians). Root Middle Leaf Total Weighted 6901.51135.11655.29691.83183.8 still more similar to MT computations, instead of a simple sum of different-level costs, we used a weighted cost c=cleaf +cmiddle 3+croot 6,(8) expressing the fact that in real problems the leaf computations are the most expensive ones, whereas the root calls are the cheapest. The denominators in (8) are taken from MT simulations, cf. see Section 5. As the comparative methods operated at the top level of accuracy, i.e., the leaf level, based on Table 5 we set the multi-start starter pool size to 50 and the SEA initial population to 70. As in the Ackley 10D case, we set the SEA genetic parameters accordingly to the corresponding HMS leaf-level parameters. We stopped the SEA right after exceeding 4000 calls. The main results of the test are shown in Table 6. The HMS successfully found all Table 6. Statistics of finding all global minima (three Gaussians). HMS Multi-start SEA Fully successful runs 30/30 16/30 0/30 Avg. number of mins 32.47 0 three minima in every execution. The multi-start failed in almost half of the runs. Moreover, in two runs it found only one global minimum. The SEA never succeeded, obtaining an average computed global minimum of 0.57 and a best found solution with the value of 0.18.As the results show, the three Gaussians benchmark is quite An agent-oriented hierarchic strategy for solving inverse problems 491 a difficult optimization test, despite its relatively low dimensionality. The source of the difficulty is the small volume of global minima’s attraction basins and the presence of significant plateaus. 4. Magnetotelluric inverse problem The solar wind produces a radiation pressure that causes a compression on the day-side and a tail on the night-side onto the magnetosphere. Due to this interaction, hydromagnetic waves are created. When those waves reach the ionosphere, they induce an EM field that works as a power source in magnetotellurics. Depending upon the type of source we are dealing with, the geomagnetic fluctuations range between the frequencies of 10−3− 105Hz, which allows us to make measurements with a resolution that ranges from a few meters to hundreds of kilometers (Vozoff, 1972) The magnetotelluric (MT) technique is used to determine a resistivity map of the Earth’s subsurface by performing electromagnetic (EM) measurements. The main difference of MT with respect to other geophysical measurement acquisition scenarios (e.g., marine controlled electromagnetic measurements) is that MT uses natural electromagnetic radiation sources generated within the ionosphere, instead of human powered antennas. Thus, acquisition of MT measurements is rather inexpensive, and can cover large areas. Applications of MT measurements include hydrocarbon (oil and gas) exploration and finding suitable regions for storage of CO2. 4.1. Forward problem. MT measurements are governed by electromagnetic phenomena, which can be described by Maxwell’s equations. When the electrical field depends only upon two spatial variables (x, z), then two independent and uncoupled modes are derived from these equations, namely, transverse electric (TE) and transverse magnetic (TM). The TE mode involves (Ey,H x,H z)field components while TM uses (Hy,E x,E z). We focus on the TE mode and we solve the equation for Ey. We decouple Maxwell’s equations by pre-multiplying both sides of Faraday’s law by μ−1 and applying the curl operator. Then, after incorporating Ampere’s law, we substitute component by component the result into the double curl operator on the electric field to obtain the equation for the y-component of the electric field: −Δu−k2u=fin Ω⊂R2,(9) where Ωis a simply connected bounded domain with a Lipschitz boundary and u(x, z)=Ey(x, z),(x, z)∈Ω. σ,(σ>σ 0>0in Ω) is the electrical conductivity field, k2=ω2μ −jμωσ,whereω=0is the wave frequency, , μ ∈R,,μ > 0stand for the permittivity and permeability of the medium considered. f=−jωμJimp y, where Jimp yis a prescribed, impressed electric current density radiating in the y-direction. To obtain the corresponding variational formulation, we pre-multiply Eqn. (9) by a test function v∈ H1 0(Ω; C). After integrating by parts and incorporating the Dirichlet boundary conditions over ΓD=∂Ω,the following abstract variational formulation (suitable for finite element computations) is obtained: Find u(σ)∈H1 0(Ω; C)such that b(σ;u(σ),v)=F(v),∀v∈H1 0(Ω; C),(10) with b(σ;u, v)=Ω ∇v∇u−Ω k2vu (11) F(v)=Ω vf, (12) where now we assume that σ∈L∞(Ω) and Jimp y∈ L2(Ω; C). The exact solution u(σ)=Eyto the problem (10) is called the primal forward solution. Let us define now the functionals Li:H1 0(Ω; C)→ Cassociated with the receiving antennas occupying Ωi domains, respectively, i=1,...,M, Li(v)= 1 meas(Ωi)Ωi v. (13) Now, we are able to define the set of dual forward problems, Find G(σ)∈H1 0(Ω; C)such that b(σ;v,G(σ)) = Li(v),∀v∈H1 0(Ω; C),(14) for each receiving antenna, i=1,...,M. 4.2. hp-FEM approximation. Using an hp-FEM finite dimensional internal approximation Vhp ⊂ H1 0(Ω; C)(see Demkowicz, 2006), we obtain the discrete versions of both primal and dual forward problems (10), (14) Find uh,p(σ)∈Vhp such that b(σ;uh,p(σ),v)=F(v),∀v∈Vhp,(15) Find Gi h,p(σ)∈Vhp such that b(σ;v,Gi h,p(σ)) = Li(v),∀v∈Vhp,(16) for each receiving antenna, i=1,...,M. To control the discretization error by performing grid refinements, we use the two-dimensional (2D) hp-adaptive algorithm described by Demkowicz (2006). It incorporates two basic components: 498 M. Smołka et al. tions), adaptive finite-element and discontinuous Petrov–Galerkin methods, multigrid solvers, image restoration algorithms, and multiphysics and inverse problems. Julen ´ Alvarez-Aramberri completed his degree in physics at the University of the Basque Country in 2007. After that, he studied for an M.Sc. in quantitative finance and an M.Sc. in mathematical modeling, statistics and computation at the University of the Basque Country. Since 2011 he has been a Ph.D. student within the Department of Applied Mathematics, Statistics and Operational Research at the University of the Basque Country under the supervision of David Pardo and within the Team Project MAGIQUE-3D at the University of Pau under the supervision of H´el`ene Barucq. He currently is in the last year of his Ph.D. working at the Basque Center for Applied Mathematics (BCAM). His Ph.D. is expected to be finished in 2015. His research interest include efficient implementation of numerical schemes, methods, and tools. In particular, he is focused on solving direct and inverse problems arising in applications of the magnetotelluric technique, which is used to retrieve information about the resistivity distribution of the Earth’s subsurface. Received: 5 September 2014 Revised: 8 November 2014