Clara Schmitt

Sébastien Rey-Coyrehourcq

Romain Reuillon

Denise Pumain

Half a billion simulations: Evolutionary algorithms and distributed computing for calibrating the SimpopLocal geographical model

Keywords: Simulation model, multi-agent system, calibration, evolutionary algorithm, geographical modelling, high performance computing, model validation

Multi-agent geographical models integrate very large numbers of spatial interactions. In order to validate those models large amount of computing is necessary for their simulation and calibration.

Here a new data processing chain including an automated calibration procedure is experimented on a computational grid using evolutionary algorithms. This is applied for the first time to a geographical model designed to simulate the evolution of an early urban settlement system. The method enables us to reduce the computing time and provides robust results. Using this method, we identify several parameter settings that minimise three objective functions that quantify how closely the model results match a reference pattern. As the values of each parameter in different settings are very close, this estimation considerably reduces the initial possible domain of variation of the parameters. The model is thus a useful tool for further multiple applications on empirical historical situations.

The challenge of calibrating geographical multi-agent models

Geographical simulation models of systems of cities are based on the assumption that the microgeographical interactions are likely to generate the emergence of stylized dynamics on the macrogeographical scale, which constitute one of the recurring universal characteristics of these complex systems [START_REF] Pumain | Theoretical principles in inter-urban simulation models: a comparison" Environment and Planning A[END_REF]. The heterogeneity and the large number of possible scenarii of implementation to represent these processes, the individual-based description of dynamics, the non-linearity of interactions, the diversity of forms of spatial relationships and the importance of the historical context are so many reasons for geographers to regularly use agentbased modelling as a support for reflection and experimentation [START_REF] Batty | Fifty Years of Urban modelling : Macro-Statics to Macro-Dynamics[END_REF][START_REF] Crooks | Key challenges in agent-based modelling for geo-spatial simulation[END_REF]Heppenstall et al, 2011;[START_REF] Sanders | Objets géographiques et simulation agent, entre thématique et méthodologie[END_REF]. We introduce here a simplified model called SimpopLocal [START_REF]Simprocess, a integrated method to validate Agent-based models in geography PhD in progress[END_REF][START_REF] Schmitt | Modeling settlement systems dynamics : From SimpopLocal to SimpopNet PhD[END_REF] that extends the Simpop family of models [START_REF] Pumain | The evolution of City systems, between History and Dynamics[END_REF][START_REF] Stonedahl | Genetic Algorithms for the Exploration of Parameter Spaces in Agent-Based Models[END_REF]. It describes the emergence of a system of settlements where the process of innovation that generates the growth dynamics of the settlements is made endogenous.

What is at stake in a successful agent-based model in geography is to assemble a set of individualbased mechanisms adapted to the level of resolution of the problem and to evaluate its ability to answer this problem with a high level of confidence [START_REF] Sargent | Verification and validation of simulation models[END_REF]. This evaluation usually requires a calibration stage [START_REF] Balci | Verification, validation, and testing[END_REF] during which the ability of the model to reproduce a specific dynamic or structure is assessed. This calibration phase is generally conducted by a trial and error method. When applied on SIMPOP2 it resulted in being fairly satisfied with the calibration obtained thanks to a hundred simulations [START_REF] Bretagnolle | Simulating Urban Networks through Multiscalar Space-Time Dynamics (Europe and the United States, 17th-20th Centuries)[END_REF]. It was thus not possible to determine if the corresponding values found for the parameters, although tedious to estimate through this manual calibration procedure (each new increment in the parameter values can disturb the dynamics), were the best possible estimation or simply denoted the existence of a local minimum in the phase space of the dynamic model. It is thus crucial to automate the process of calibration and to substitute temporarily the human expertise of the model by an automated expertise, which requires a quantitative transposition of the expert evaluation.

To carry out this automation, we developed and implemented methods and tools within the community of research on the complex systems [START_REF] Bourgine | French Roadmap for complex Systems 2008-2009[END_REF] and more precisely in the modelling platform called SimProcess [START_REF] Reuillon | OpenMOLE, a workflow engine specifically tailored for the distributed exploration of simulation models[END_REF]. After having described the geographical model Simpoplocal on which the method is tested, we present an automated procedure of calibration that allows a massive exploration of the value space of the parameters, which was up to now impossible to handle in the usual way. This exploration is guided by specific and relevant objectives as for the issues dealt with by SimpopLocal. After extraction and study of the results, we identify a set of coherent parameter settings that give particularly interesting simulation outputs in compliance with the model calibration objectives. The main contribution of the paper is to provide a re-useable method for calibrating a real-world ABM model.

The SimpopLocal model

A model to simulate the emergence of structured urban settlement system

SimpopLocal is a stylized model describing an agrarian society in the Neolithic period, during the primary "urban transition" manifested by the appearance of the first cities [START_REF] Schmitt | Modeling settlement systems dynamics : From SimpopLocal to SimpopNet PhD[END_REF]. It is designed to study the emergence of a structured and hierarchical urban settlement system by simulating the growth dynamics of a system of settlements whose development remains hampered by strong environmental constraints. This exploratory model seeks to reproduce a particular structure of the Rank-Size distribution of settlements well defined in the literature as a generalized stylized-fact: for any given settlement system throughout the time and continents, the distribution of sizes is strongly differentiated, exhibiting a very large number of small settlements and a much smaller number of large settlements [START_REF] Archaeomedes | Des oppida aux métropoles[END_REF][START_REF] Berry | Cities as systems within systems of cities[END_REF][START_REF] Fletcher | Settlement Archaeology: World-Wide Comparisons[END_REF][START_REF] Liu | Settlement patterns, chiefdom variability, and the development of early states in North China[END_REF]. This distribution can be modeled by a power-law when small settlements are not considered or a lognormal when small settlements are considered [START_REF] Favaro | Gibrat Revisited: An Urban Growth Model including Spatial Interaction and Innovation Cycles[END_REF][START_REF] Robson | Urban growth, an approach[END_REF]. Such distributions are easily simulated by rather simple and non spatial statistical stochastic models [START_REF] Gibrat | Les inégalités économiques[END_REF][START_REF] Simon | On a class of Skewed Distribution Functions[END_REF]. But for theoretical considerations and to specify the model in a variety of geographical frames, we think it necessary to make explicit the spatial interaction processes that generate the evolution of city systems. According to the evolutionary theory of system of cities [START_REF] Pumain | The evolution of City systems, between History and Dynamics[END_REF], the growth dynamics of each settlement is controlled by its ability to generate interurban interactions. The multi-agent system modeling framework enables us to include mechanisms, derived from this theory, that govern the interactions between settlements [START_REF] Bretagnolle | From theory to modelling : urban systems as complex systems[END_REF][START_REF] Pumain | Theoretical principles in inter-urban simulation models: a comparison" Environment and Planning A[END_REF]. The application of this concept resulted in several Simpop models [START_REF] Bretagnolle | Simulating Urban Networks through Multiscalar Space-Time Dynamics (Europe and the United States, 17th-20th Centuries)[END_REF]Bura et al, 1996;[START_REF] Pumain | The Future of Urban Systems[END_REF][START_REF] Sanders | SIMPOP: a multi-agent system for the study of urbanism[END_REF] in which the expected macro-structure of the log-normal distribution of sizes emerges from the differentiated settlement growth dynamics induced by the heterogeneous ability of interurban interactions. Therefore the aim of SimpopLocal is to simulate the hierarchy process via the explicit modeling of a growth distribution that is not entirely stochastic as in the Gibrat model [START_REF] Gibrat | Les inégalités économiques[END_REF] but that emerges from the spatial interactions between micro-level entities. Compared to the previous Simpop models, the originality of SimpopLocal is twofold. First, it is a simplified version which no longer qualitatively distinguishes the successive urban functions simulated during the evolution of urban systems but transposes them into an abstract innovation process having over time less and less impact in terms of wealth and population growth (this rule is in accordance with the concept of the decreasing efficiency of the improvement of an existing socio-economic system and prevents introducing any a priori amplification effect). Second, the appearance process of these innovations is an endogenous process that is linked to the size of the settlements.

Agents, attributes, and mechanisms of the model

The SimpopLocal model is a multi-agent model developed with the Scala programming language. It simulates the growth dynamics of agrarian settlements and their possible evolution toward urban settlements under strong environmental constraints that are progressively overcome by successive innovations.

The landscape of the simulation space is composed of hundred of settlements. Each settlement is considered as a fixed agent and is described by three attributes: the location of its permanent habitat, the size of its population, and the available resources in its local environment. The amount of available resources is quantified in units of inhabitants and can be understood as the carrying capacity of the local environment for sustaining a population which depends on the resources exploitation skills that the local population has acquired from inventing or acquiring innovation. Each new innovation acquired by a settlement develops its exploitation skills. This resource exploitation is done locally and sharing or trade is not explicitly represented in the model.

The growth dynamics of a settlement is modelled according to the assumption that its size is dependent on the amount of available resources in the local environment and is inspired from the Verhulst model [START_REF] Verhulst | Recherches mathématiques sur la loi d'accroissement de la population[END_REF] or logistic growth. For this experiment, we assume a continuous general growth trend for population -this may be different in another application of the model. The growth factor r is expressed on an annual basis, thus one iteration or step of the model simulates one year of demographic growth. The limiting factor of growth is the amount of available resource that depends on the number M of innovations the settlement i has acquired by the end of the simulation step t.

is the population of the settlement i at the time t :

= 1 + (1 - )
The acquisition of a new innovation by a settlement allows it to overtake its previous growth limitation by allowing a more efficient extraction of resources and thus a gain in population size sustainability. With the acquisition of innovations the amount of available resources tends to the maximal carrying capacity R max of the simulation environment:

The mechanism of this impact follows the Ricardo model of diminishing returns (also a logistic model (Turchin, 2003)). The InnovationImpact represents the impact of the acquisition of an innovation and has a decreasing effect on the amount of available resources with the acquisition of innovations:

= 1 + !"#$%! &'#($(1 - )
Acquisition of innovations can occur in two ways, either by the emergence of innovation within a settlement or by its diffusion through the settlement system. In both cases, interaction between people inside a settlement or between settlements is the driving force of the dynamics of the settlement system. It is a probabilistic mechanism, depending on the size of the settlement. Indeed, innovation scales super-linearly: the greater the number of innovations acquired, the bigger the settlement and the higher the probability of innovating [START_REF] Arthur | The Nature of Technology: What It Is and How It Evolves[END_REF]Diamond, 1999;Lane et al, 2009). To model the super-linearity of the emergence of innovation within a settlement, we model its probability by a binomial law. If P creation is the probability that the interaction between two individuals of a same settlement is fruitful, i.e. leads to the creation of an innovation, and N the number of possible interactions, then by the Binomial law, the probability of the emergence of at least one innovation P(m creation >0) can be calculated and then used in a random drawing :

(& )* > 0) = 1 -(& )* = 0) = 1 -- .! 0!(.10)! )* 0 (1 -)* ) .10 2 = 1 -(1 -)* ) .
If the size of the settlement is then the number N of possible interactions between individuals of that settlement is :

3 = ( -1) 2
The diffusion of an innovation between two settlements depends both on the size of populations and the distance between them. If P diffusion is the probability that the interaction of two individuals of two different settlements is fruitful, i.e. leads to the transmission of the innovation, and K the number of possible interactions, then by the Binomial law the probability of diffusion of at least one innovation P(m diffusion >0) can be calculated and used in a random drawing :

(& 5 66 > 0) = 1 -(1 -5 66 ) 7
But in this case, the size K of the total population interacting is a fraction of the population of the two settlements i and j which is decreasing by a factor DistanceDecay with the distance D ij between the settlements, as in the gravity model [START_REF] Wilson | A family of spatial interaction models, and associated developments[END_REF]):

8 = 9 2 : 9 ; *;* <
The process of population growth and the process of innovation creation/diffusion are re-iterated throughout the simulation. Because of the two positive feedbacks that operate on resource and population growth through the creation of innovation, the model is able to generate a very rapid expansion of settlements, an escalation of settlement growth. The simplest way to avoid situations where innovations are created in too large numbers which would lead to huge time consuming simulations, is to decide to stop the simulation when it reaches an arbitrary number of 10 000 innovations. In order to ensure the replicability of the model, the source code of SimpopLocal is filed in a public repository: http://iscpif.github.io/simpoplocal-epb/.

Parameters to calibrate

SimpopLocal has many parameters that have to be estimated in order to calibrate the model. Some can be empirically evaluated with the help of historical data and knowledge, while it is very difficult to give values to others. Those regarding the initial spatial distribution and organisation of settlements in the landscape can be approximated. The log-normal distribution of the settlement sizes and the central place theory [START_REF] Christaller | Die Zentralen Orte in Süddeutschland. Eine ökonomisch-geographische Untersuchung über die Gesetzmäβigkeit der Verbreitung und Entwicklung der Siedlungen mit städtischen Funktionen[END_REF] for the geographical distribution of locations are models that are widely used by archaeologists to describe their spatial data [START_REF] Archaeomedes | Des oppida aux métropoles[END_REF][START_REF] Johnson | Aspects of regional analysis in archaeology[END_REF][START_REF] Sanders | Regards croisés de géographes, économistes et archéologues sur la hiérarchie des systèmes de peuplement : de l'empirie aux systèmes complexes[END_REF], including Neolithic archaeological sites [START_REF] Liu | Settlement patterns, chiefdom variability, and the development of early states in North China[END_REF]. In SimpopLocal, the mean density of that landscape and the average size of each settlement are representative of the usual orders of magnitude presented in these works. A hundred settlements are distributed according to these two theories and each settlement is initially composed of some 80 up to 400 inhabitants. Several scholars agree that an average annual growth of 0.02% is representative of the growth of agrarian settlements in the Neolithic times [START_REF] Bairoch | De Jéricho à Mexico: villes et économie dans l'histoire[END_REF][START_REF] Pumain | Theoretical principles in inter-urban simulation models: a comparison" Environment and Planning A[END_REF]. The length of time required for a transition from agrarian to urban settlements is estimated according to [START_REF] Bairoch | De Jéricho à Mexico: villes et économie dans l'histoire[END_REF]Marcus and Sabloff, 2008) to a couple of thousand years. We choose to operate our simulations on a four thousand year time period for settlements ranging from one hundred inhabitants up to about ten thousand inhabitants.

Because of a lack of empirical data, five parameters cannot be empirically approximated:

-P creation , the probability that an innovation emerges from the interaction between two individuals of a same settlement.

-P diffusion , the probability that an innovation is transmitted between two individuals of different settlements. We consider that the probability of diffusion is greater than the probability of creation, which means that copying is easier than inventing [START_REF] Pennisi | Conquering by Copying[END_REF] in the model.

-

InnovationImpact, the impact of the acquisition of innovation on the growth of settlements.

-DistanceDecay, the deterrent effect of distance on diffusion.

-R max , the maximum carrying capacity of the landscape of each settlement (measured in number of inhabitants).

The model has been simplified as far as possible to retain only five parameters that cannot be empirically estimated. This number may seem rather low compared to the 40 or 50 parameters that were activated in the other versions of SIMPOP models, but this simplification allows a global exploration of the capabilities of SimpopLocal. Indeed, by means of intensive exploration, a calibrated state of the model can produce an estimation of the value of some parameters that could not be deduced from the empirical literature on archaeological settlements. These estimations will be useful afterwards for making predictions about the possible evolution of concrete settlements systems or comparing the evolutions of early urban systems in different regions. SimpopLocal is conceived for providing an open evolution that may lead to any type of size distribution of settlements depending on the parameter values. In order to estimate parameter values that could generate plausible settlement size distributions, we designed an automated calibration procedure.

Designing an automatised calibration procedure

Calibration as an optimisation problem

Model calibration is a procedure which seeks to minimize the difference between the behaviour simulated by the model and a behaviour defined according to expert knowledge and/or data. In most multi-agent systems the calibration is done manually by introducing values for the parameters and visually verifying that the output of the model corresponds to the expected results. But this method raises numerous problems highlighted by Stonedahl [START_REF] Stonedahl | Genetic Algorithms for the Exploration of Parameter Spaces in Agent-Based Models[END_REF]: relationships between parameters are often non-linear; expertise during behaviour exploration may be biased by erroneous assumption; some parameters cannot be compared to empirical values; and manual exploration is tedious. In our first attempt to manually calibrate the model, we could not reach a stage where the calibration could reasonably be considered satisfying.

Because of the number of parameters and of the continuous scale of variation of their fields of variation, any exhaustive search strategy is not tractable either. First of all, this strategy generates a number of experiments that grows exponentially with the number of parameters under study. This would imply too gigantic an amount of computation in the case of SimpopLocal. Moreover, that strategy also produces a large quantity of data that then has to be processed and visualized. The search for patterns in the space of the dynamics of the model is thus replaced by the problem of a search for patterns in a database. Post-processing of such a large quantity of data is long and tiresome. Instead of resorting to exhaustive search methods, we choose to consider calibration as an optimisation exercise. The a posteriori exploration of results becomes then an a priori question: is there a parameter setting that would match our expectation?

Evolutionary algorithms have been established as a suitable solution to this way of considering calibration [START_REF] Calvez | Automatic tuning of agent-based models using genetic algorithms[END_REF][START_REF] Stonedahl | Genetic Algorithms for the Exploration of Parameter Spaces in Agent-Based Models[END_REF]. They have been used to calibrate multi-agent systems in several application fields such as medicine [START_REF] Castiglione | Optimization of HAART with genetic algorithms and agent-based models of HIV infection[END_REF], ecology [START_REF] Duboz | Application of an evolutionary algorithm to the inverse parameter estimation of an individual-based model[END_REF], economics [START_REF] Espinosa | A Genetic Algorithm for the Calibration of a Micro-Simulation Model[END_REF]Stonedahl et al., 2010(a)), hydrology [START_REF] Solomatine | Automatic calibration of groundwater models using global optimization techniques[END_REF]... Despite the wide use of multi-agent systems in social sciences, this method has not yet been applied very often. To our knowledge only a few real world applications have been carried out [START_REF] Heppenstall | Genetic algorithm optimisation of an agent-based model for simulating a retail market[END_REF]Stonedahl et al., 2010(b)). Indeed this kind of numerical experiment remains a real challenge:

-it requires the definition of quantitative goals that evaluate the simulation outputs in a coherent manner which correspond to experts' expectations, -it generates a massive computation load that requires technical skills and adapted infrastructures, -it seeks to optimise noisy (stochastic) fitness function which is in itself a challenging exercise [START_REF] While | Applying evolutionary algorithms to problems with noisy, time-consuming fitness functions[END_REF].

To overcome these obstacles, cutting edge knowledge and tools in several fields of expertise have to be coupled within highly transdisciplinary teams involving social scientists, statisticians, optimization method specialists and distributed-computing experts.

Exploring the parameter space given a set of objectives

Evolutionary algorithms are heuristics that scan the search space using strategies inspired by natural processes to solve optimisation problems. In order to use evolutionary algorithms to calibrate a model, the first step is to formalize what the expected result of a "good" simulation is, which is, in our case, a suitable configuration of the settlement system. This suitable configuration reflects stylized facts that were established over the years by many scholars thanks to the exploration and processing of large amounts of empirical data. We have summarized this suitable configuration by three objective functions extracted as relevant stylised facts from the existing literature: a log-normal distribution of settlement sizes [START_REF] Archaeomedes | Des oppida aux métropoles[END_REF][START_REF] Johnson | Aspects of regional analysis in archaeology[END_REF][START_REF] Liu | Settlement patterns, chiefdom variability, and the development of early states in North China[END_REF][START_REF] Sanders | Regards croisés de géographes, économistes et archéologues sur la hiérarchie des systèmes de peuplement : de l'empirie aux systèmes complexes[END_REF], a maximum size of ten thousands inhabitants [START_REF] Bairoch | De Jéricho à Mexico: villes et économie dans l'histoire[END_REF]Marcus and Sabloff, 2008), and a four thousand years period duration for achieving this distribution (in any region of the world where agriculture was invented, the emergence of the first cities occurred a few thousand years later [START_REF] Bairoch | De Jéricho à Mexico: villes et économie dans l'histoire[END_REF]Marcus and Sabloff, 2008).

The automatic technique of calibration that we propose is based on an evolutionary algorithm that explores the space of the parameter settings. This exploration is guided by these three objectives [START_REF]Simprocess, a integrated method to validate Agent-based models in geography PhD in progress[END_REF][START_REF] Schmitt | Modeling settlement systems dynamics : From SimpopLocal to SimpopNet PhD[END_REF]. Each parameter setting (i.e. set of values for the 5 unknown parameters of the model) is evaluated according to the simulation outputs it produces. This evaluation measures the proximity between the outputs of simulation and the three objective functions defined for the model. This evaluation thus assesses the ability of the parameter settings to reproduce the stylized facts that we seek to simulate. The parameter settings receiving the best evaluations are then used as a basis for generating new parameter settings which will then be tested.

The SimpopLocal model being stochastic, the simulation outputs do vary for a given parameter setting from one simulation to the next. Therefore, the evaluation of the parameter setting on the three objectives must take into account this variability. We checked if 100 simulations for each set of parameters was sufficient to capture this variability and found out that it was a suitable compromise between capturing the variability and not increasing too much the computation duration. Here is how we test each set of parameters according to the three objective functions:

-The objective of distribution quantifies the ability of the model to produce settlement size distributions that fit a log-normal distribution. First we evaluate the outcome of each subset of 100 simulations corresponding to one parameter setting by computing, according to a 2-sample Kolmogorov-Smirnov test, the deviation between the simulated distribution and a theoretical lognormal distribution having the same mean and standard deviation. Two criteria are reported, with value 1 if the test is rejected and 0 otherwise: the likelihood of the distribution (the test returns 0 if p-value > 5%) and the distance between the two distributions (the test returns 0 if D-value < D α 1

). In order to summarize those tests in a single quantified evaluation, we add the results of the two tests on the 100 simulations. The best possible score on this objective is thus 0 (all tests returning 0) and the worst 200 (all tests returning 1). By dividing this score by the worst possible value (200), we get a normalised error that facilitates the comparisons between parameter settings according to their scores on the three objectives.

-The objective of population quantifies the ability of the model to generate large settlements that have an expected size. The outcome of one simulation is tested by computing the deviation between the size of the largest settlement and the expected value of 10000 inhabitants: |(population of largest settlement -10000)/10000|. The evaluation of the parameter setting reports the median of this test on the 100 simulations. This value represents the normalised error produced by the parameter setting being evaluated on the calibration (the closest to 0, the smallest the error). Computing this normalised error is not necessary but it facilitates the comparisons between parameter settings according to their scores on the three objectives.

-The objective of simulation duration quantifies the ability of the model to generate an expected configuration in a suitable length of time (in simulation steps). The duration of one simulation is tested by computing the deviation between the number of iterations of the simulation and the expected value of 4000 simulation steps: |(simulation duration -4000)/4000|. The evaluation of the parameter setting reports the median of this test on the 100 simulations. This value represents the normalised error produced by the parameter setting being evaluated on the calibration. This normalised error (the closest to 0, the smallest the error) facilitates the comparisons between parameter settings according to their scores on the three objectives.

The three objective functions constitute a fitness function of our optimisation problem which is therefore multi-objective. To solve such problems, it is sometimes possible to aggregate the vector values in a single scalar value reflecting an absolute quality of the solution. In our case it is not possible to provide a meaningful aggregated function. Thus we rely on a multi-objective optimisation algorithm. This type of algorithm computes compromise solutions: among those solutions none dominates the others by presenting better scores on all three objectives at the same time. Such solutions are called Pareto compromise and the entire set of Pareto compromise solutions is called the Pareto front.

Executing a large scale distributed evolutionary algorithm with a noisy time-consuming fitness function

The use of global methods of search like evolutionary algorithms for the calibration of a multi-agent model (and especially a stochastic multi-agent model) entails a high computational cost [START_REF] Sharma | Auto-calibration of hydrological models using high performance computing[END_REF]. This kind of load is too large to be executed on local computers. Super-computers are very expensive and not easily available in most laboratories, therefore we decided to develop a calibration framework that could run on both expensive super-computers and on computing grids which make computing power more accessible by federating it on a world wide scale. For the experiment presented in this paper we used the European Grid Infrastructure (EGI).

Computing grids offer a solution for the resolution of such computationally intensive problems, however computing at such a large scale is a true challenge. It supposes orchestrating the execution of tens of thousands of instances of the model on computers distributed all over the world. The cumulated probability of local breakdowns and the impossibility to optimally distribute the workload on the grid system makes its effective use very difficult. To overcome these difficulties we have used the OpenMOLE platform (Open Model Exploration, www.openmole.org) [START_REF] Reuillon | Declarative task delegation in OpenMOLE[END_REF][START_REF] Reuillon | OpenMOLE, a workflow engine specifically tailored for the distributed exploration of simulation models[END_REF], which provides a bridge that helps modellers to cross the technical and methodological gap which separates them from high performance computing. OpenMOLE is a dedicated, textual and graphic language, exposing coherent bricks at the right level of abstraction to conceive reusable experiments in order to solve inverse problems on models. The methods of resolution are described independently of a particular model and thus support the reproducibility and the re-use of the experimental numerical methods between the modellers.

For calibrating SimpopLocal we implemented a classical multi-objective evolutionary algorithm called SMS-EMOA (Emmerich 2005) in a distributed manner, on top of OpenMOLE. The parameters of this algorithm are exposed in table 1. The code of this workflow and the code for the fitness computation are available on: http://iscpif.github.io/simpoplocal-epb/. This algorithm evaluates millions of parameter settings to approximate the Pareto front. In its most efficient implementation one execution of SimpopLocal lasts 1.5 seconds on a current processor. The automatic calibration of SimpopLocal requires the simulation of almost 500 million executions of the model, which would represent nearly 20 consecutive years on a single computer. To achieve this huge computation load, the SMS-EMOA algorithm was distributed on the numerous computers of the EGI by using the technique known as the "island model" (Emmerich 2005;[START_REF] Whitley | Island model genetic algorithms and linearly separable problems[END_REF].

The classical island models consist in instantiating permanent islands (isolated instances of an evolutionary algorithm) on many computers and organising migration of solutions between those islands. The EGI grid is a worldwide batch system on which organising direct communications between islands running on multiple execution nodes is very challenging. Thus the classical island model has been adapted to the EGI architecture.

During the computation, a central population of 200 parameter settings (here called individuals) is maintained on the computer that orchestrates the submission of the computing jobs on the grid. Each job computes the evolution of the population of an island, which is an independent instance of SMS-EMOA seeded with 50 individuals randomly sampled among the central population of 200 individuals. The "island job" life cycle is managed by the EGI : each job is submitted to the EGI, then to the queue of a cluster and finally starts running when a slot becomes available on the cluster.

When it starts running it is configured to run for 1 hour. 5000 concurrent jobs are maintained on the grid at any time.

Once an island job is finished its final population of individuals (i.e. parameter settings) is transferred back to the submission computer and merged into the global population using the elitism algorithm of the SMS-EMOA (based on the contribution to the hypervolume of the Pareto). A new island job is then submitted to the grid. This algorithm is run until 190 000 island jobs have been completed.

As previously stated, the fitness function of this experiment is stochastic, so we execute the model a hundred times for each parameter setting evaluation. Despite this precaution, it is well-known that among the millions of executions computed by the algorithms, some solutions might get "lucky" and be over-estimated, creating a bias in the evolution process [START_REF] While | Applying evolutionary algorithms to problems with noisy, time-consuming fitness functions[END_REF]. According to Pietro et al. ( 2004) who reviewed various method to tackle this problem, we choose to re-evaluate already evaluated solutions on a regular basis to prevent "from retaining incorrect fitness" (p.2). A specific strategy was implemented: one over a hundred parameter setting evaluations was used to reevaluate already evaluated parameter settings instead of evaluating a new parameter setting. A parameter setting which has been over-estimated the first time and thus kept among the optimal individuals by mistake, will probably be correctly evaluated the second time and eliminated from the optimal selection.

Results and discussion

Evolutionary algorithms are good solutions when no exhaustive method is available, however they don't ensure the optimality of the computed solution. Because of their heuristic nature, they are procedures that can run indefinitely. In this work we have used the number of evaluated islands as a stopping criterion given the computing power that was at our disposal. That way we preconditioned the experience so that the total execution would last for 3 to 4 days, which means that the evolutionary optimisation algorithm was stopped after the execution of 190000 islands, which represents about 200000 computation hours (200000 CPU hours 2 ). An evolutionary algorithm is declared "converged" when it makes no further improvements in the search for good solutions. One of the best metrics presently available to measure the convergence of multi-objective optimisation algorithm is the stagnation of the hypervolume [START_REF] Fonseca | An improved dimension-sweep algorithm for the hypervolume indicator[END_REF][START_REF] Naujoks | Multi-objective optimisation using S-metric selection: Application to three-dimensional solution spaces[END_REF][START_REF] Zitzler | Multiobjective optimization using evolutionary algorithms-A comparative case study[END_REF]. The hypervolume measures the volume of the dominated portion of the objective space and its stagnation indicates that the algorithm has converged. Figure 1 depicts the evolution of the hypervolume in function of the number of evaluated islands, which can also be quantified in terms of computation time or CPU hours. After the execution of 120000 CPU hours, the hypervolume stabilises. This stability in the long term (70000 CPU hours) indicates that the algorithm has converged. The Pareto front probably wouldn't have improved much more with further computation time. This result shows that despite the stopping criterion that we have chosen because of technical limitations, the proposed parameter settings of the Pareto front are seemingly the best possible results that can be obtained with this method and they most likely correspond to the global optimum.

Figure 1 : Hypervolume enclosed by the Pareto front and the converging evolution of the evolutionary algorithm under the constraint of 3

After the convergence of the algorithm, calibrated sets of parameter value approximation of the objective scores they are well evaluated, 10000 conducted. After re-computing the dominance selection, excluded and 62 parameter settings exploration does not lead to the selection of a single optimised solution candidates for the calibration: each set of specific compromise on the three objectives. settings those satisfying only compromise on the three objectives objective. This configuration confirms that these procedure vacillates between them to find well settings of the Pareto front, some are less satisfying than others. parameter setting that have at least error over 0.1, i.e. over 10% of error values estimated for each parameter enclosed by the Pareto front and the converging evolution of the algorithm under the constraint of 3 objectives 3 . convergence of the algorithm, 200 different parameter settings were proposed as possib sets of parameter values. Because a set of 100 executions can only provide a approximation of the objective scores of a parameter setting and because we want to 10000 more executions of each proposed parameter settings where computing the dominance selection, the dominated parameter settings parameter settings remained in the Pareto front. As said above to the selection of a single optimised solution but leads to a set of possible candidates for the calibration: each set of parameter settings selected by the procedure represents a compromise on the three objectives. It is thus possible to distinguish among the parameter two objectives out of three from those which offer a better jectives without however reaching the best possible values confirms that these three objectives are not redundant them to find well-fitted parameter settings. Among the , some are less satisfying than others. We consider as uns at least one of the scores of the three objectives that of error. Only 29 parameter settings satisfy this new condition. values estimated for each parameter of the 29 selected sets are presented in table The analysis of this new subset shows interesting aspects. First, all of them correctly identify the order of magnitude of the only parameter of the model which can be roughly deduced from the initial conditions and the calibration objectives: the maximum resource parameter (R max ) which is used to limit the logistic growth of the population of each settlement5 . Second, the four parameters whose values are a priori unpredictable (DistanceDecay, the parameter of dissuasion from the interactions by the distance; P creation , the probability of appearance of an innovation; P diffusion , the probability of its diffusion; and finally InnovationImpact, the impact of the innovation on the growth of the population) are estimated in a very small domain of variation compared with their possible domain of variation (table 3). The ratio between the estimated and theoretical volume defined by these five dimension domains is of ~ 4.7 -16 / 320000 = 1.5 -21 ! What is remarkable is how the values taken by each parameter are all in a small neighbourhood, which suggests high reliability considering that they must be given those orders of magnitude to obtain plausible results while using the model for simulation. However it should be noted that the exact values estimated for each parameter do not have any absolute meaning. They only make sense all together, according to their interrelationships in the mechanisms of the model. Within the 29 parameter settings subset there is no way to prefer a setting over the others: the multiple executions of each setting lead to good results and acceptable variability of the outputs. Figure 2 and3 give an example of such outputs. They depict the results produced by one of the parameter settings (in bold characters6 in table 2). Figure 2 shows the evolution of the rank-size distribution of settlement sizes during a simulation of the parameter setting. This evolution corresponds to what is expected of the model: a progressive and continuous process of hierarchical organisation of the settlement system (the slope of the linear fit of the rank-size distribution shifts from 0.2 to 0.9 in 4000 years for a maximum size reached of about 10000 inhabitants). This result is quite robust if we consider the low variability of the recorded final state from one simulation to another (figure 3). This paper illustrates the potential of a new modelling procedure for exploring the dynamics of geographical systems which may be useful as well for other social sciences combination of three principles: a conceptual framework making explicit and quantifying the criteria for evaluating the model; a methodology for exploring the parameter space with an automated process using evolutionary algorithms;

Rank-Size distribution at the final simulation step parameter settings.

procedure used in this article was applied for the first time agent model intended to simulate the emergence of a system of urban settlements convincing results to solve the calibration problem for multi whose value could not be estimated from empirical data, it would have been find an estimation of a calibrated parameter setting "manually". By reversing this calibration problem into an optimisation problem, this automated calibration procedure generates gs which reproduce the stylized facts embodied in three objective functions very enables us to confirm a decisive advance in the validation of this model : the SimpopLocal model, as conceived in its simplicity, is able to satisfactorily generate a plausible evolution including the emergence of an urban hierarchy.

This paper illustrates the potential of a new modelling procedure for exploring the dynamics of geographical systems which may be useful as well for other social sciences. The novelty is in the combination of three principles: a conceptual framework making explicit and quantifying the criteria for evaluating the model; a methodology for exploring the parameter space with an automated algorithms; a technical protocol for reducing the at the final simulation step on 100 simulations procedure used in this article was applied for the first time on a multiagent model intended to simulate the emergence of a system of urban settlements called esults to solve the calibration problem for multi-agents models. , it would have been parameter setting "manually". By reversing this calibration problem into an optimisation problem, this automated calibration procedure generates gs which reproduce the stylized facts embodied in three objective functions very the validation of this model : we proved that torily generate a plausible This paper illustrates the potential of a new modelling procedure for exploring the dynamics of . The novelty is in the combination of three principles: a conceptual framework making explicit and quantifying the criteria for evaluating the model; a methodology for exploring the parameter space with an automated a technical protocol for reducing the exploration time by massive distributed computing. The application of these principles to modelling not only improves the quality of communication about the model, but as well ensures the repeatability of its experimentation and contributes to create accessible modelling tools that can be shared via the free and open source software tool for model experiments, OpenMOLE (www.openmole.org). In order to help the dissemination of these techniques, we published several tutorials on the developed tools : on how to design a grid exploration of a NetLogo model 7 , how to explore a NetLogo model with evolutionary algorithms 8 , etc. This work forges a path for other modelers, who could use a similar data-intensive grid-based approach in their modeling experiments.

Figure 2 :

 2 Figure 2: Evolution of the Rank-Size distribution during a simulation of one of the best calibrated parameter settings.

Figure 3 :

 3 Figure 3: Variability of the Rank of one of the best calibrated parameter settings

Table 1 . Parameters of the implemented evolutionary algorithm

 1

	Evolutionary Algorithm	SMS-EMOA

Table 2 . Calibrated Parameter settings

 2

		PARAMETER SETTING
	R max	Distance-	P creation
		Decay
	10464	0.691 1.11E-06
	10465	0.709 1.14E-06
	10459	0.705 1.15E-06
	10465	0.708 1.15E-06
	3 In order to smooth the local variability of the raw hypervolume
	island sliding window.

4

For readability the values presented are rounded.

. Calibrated Parameter settings 4

  In order to smooth the local variability of the raw hypervolume value, it has been averaged using a 1000For readability the values presented are rounded.

	PARAMETER SETTING		% ERROR ON EACH
	P diffusion	Innovation-	%error	%error
		Impact	Distribution-	Population
			Objective	Objective
	7.90E-07	0.0077	2.02
	7.88E-07	0.0078	2.04
	7.92E-07	0.0077	2.22
	7.86E-07	0.0078	2.51

enclosed by the Pareto front and the converging evolution of the

	10261	0.679 1.17E-06	7.78E-07	0.0078	2.79	5.29	4.55
	10262	0.679 1.15E-06	7.52E-07	0.0078	2.95	5.00	2.15
	10262	0.665 1.14E-06	7.24E-07	0.0078	2.99	5.53	1.23
	10261	0.683 1.12E-06	7.38E-07	0.0080	3.1	4.34	1.15
	10260	0.699 1.17E-06	7.38E-07	0.0079	3.62	5.51	0.20
	10287	0.690 1.23E-06	7.56E-07	0.0078	3.63	3.74	3.46
	10259	0.688 1.20E-06	7.41E-07	0.0079	3.74	3.55	2.48
	10169	0.736 1.29E-06	7.39E-07	0.0079	4.90	5.20	0.03
	10205	0.683 1.19E-06	7.42E-07	0.0082	5.10	2.60	6.03
	10126	0.738 1.22E-06	7.61E-07	0.0082	5.69	2.88	1.50
	10126	0.738 1.24E-06	7.39E-07	0.0082	6.02	3.08	0.55
	10096	0.701 1.14E-06	7.14E-07	0.0084	6.12	2.58	1.55
	10169	0.736 1.29E-06	7.39E-07	0.0080	6.25	2.46	1.20
	10165	0.734 1.29E-06	7.24E-07	0.0080	6.31	e proposed as possible 2.91 0.30
	10121	0.732 1.28E-06	7.41E-07	0.0081	100 executions can only provide an 6.41 2.36 1.90
	10164	0.735 1.29E-06	7.27E-07	0.0080	and because we want to make sure that 6.45 2.74 0.45
	10103	0.733 1.24E-06	7.42E-07	0.0084	proposed parameter settings where 7.67 1.90 3.10
	10092	0.736 1.29E-06	7.14E-07	0.0082	7.81	parameter settings were 2.22 1.10
	10098	0.737 1.29E-06	7.12E-07	0.0082	7.84	above, a multi-objective 2.55 0.58
	10094 10129 10110 10091 10091 10109	0.741 1.28E-06 0.737 1.29E-06 0.735 1.28E-06 0.744 1.31E-06 0.741 1.31E-06 0.734 1.28E-06	but leads to a set of possible 8.46 1.99 1.00 selected by the procedure represents a 8.64 1.97 0.68 to distinguish among the parameter 9.04 2.48 0.03 those which offer a better 9.22 1.51 2.68 without however reaching the best possible values on each 7.12E-07 0.0083 7.07E-07 0.0082 6.77E-07 0.0083 7.25E-07 0.0083 7.12E-07 0.0083 9.61 1.49 2.15 redundant and that the mong the 62 parameter 6.79E-07 0.0084 9.64 1.77 0.18
					We consider as unsuited the
						that has a normalised
						this new condition. The
					presented in table 2.
						EACH OBJECTIVE SCORE
						%error	%error
						Population-	Time-
						Objective	Objective
						6.59	1.00
						7.31	0.10
						6.98	0.48
						6.00	0.40
					value, it has been averaged using a 1000

Table 3 . Global domain of exploration and calibrated parameter value domain of variation

 3

		Explored domain of variation	Domain enclosed in the Pareto
			front
	P creation	[0 ; 1]	~[ 1.1E-06; 1.3E-06 ]
	P diffusion	[0 ; 1]	~[ 6.7E-07; 6.9E-07 ]
	InnovationImpact	[0 ; 2]	~[ 7.7E-03; 8.4E-03 ]
	DistanceDecay	[0 ; 4]	~[ 0.66; 0.75 ]
	R max	[1 ; 40 000]	~[ 10090; 10465 ]
	Domain Volume	320000	~ 4.7E-16

Computed with α = 1.36.

Non-linear effects make this estimation tricky: the expected value would have been 10000 inhabitants, according to the objective of population, but in order to reach this value, the maximum carrying capacity of the landscape must be slightly higher.

The non-rounded parameter setting is : {R max : 10259, 331894632433; DistanceDecay : 0, 6882107473716844; P creation : 1.2022185310640896E-06; P diffusion :

7.405303653131592E-07; InnovationImpact : The seed used for the simulation depicted on figure2 is: -6863419716327549772.

Acknowledgments

This work is funded by the ERC Adv Grant Geodivercity, the Agence de l'Environnement et la Maitrise de l'Energie (ADEME), the network Réseau de Recherche sur le Développement Soutenable (R²DS), the Institut des Systèmes Complexes Paris Ile-de-France (ISC-PIF). Results obtained in this paper were computed on the biomed and the vo.complex-system.eu virtual organization of the European Grid Infrastructure (http://www.egi.eu). We thank the European Grid Infrastructure and its supporting National Grid Initiatives (France-Grilles in particular) for providing the technical support and infrastructure. We also thank an anonymous reviewer for his/her help in improving the first version of the paper.


Gérôme Bovet
email: gerome.bovet@telecom-paristech.fr

Gautier Briard
email: gautier.birard@utbm.fr

Jean Hennebert
email: jean.hennebert@hefr.ch

A Scalable Cloud Storage for Sensor Networks

Keywords: H.2.4 [Information Systems]: Distributed databases Cloud storage, Web-of-Things, Sensor networks

Data storage has become a major topic in sensor networks as large quantities of data need to be archived for future processing. In this paper, we present a cloud storage solution benefiting from the available memory on smart things becoming data nodes. In-network storage reduces the heavy traffic resulting of the transmission of all the data to an outside central sink. The system built on agents allows an autonomous management of the cloud and therefore requires no human in the loop. It also makes an intensive use of Web technologies to follow the clear trend of sensors adopting the Web-of-Things paradigm. Further, we make a performance evaluation demonstrating its suitability in building management systems.

INTRODUCTION

Since a few years, Web technologies are gaining more importance for communication between so-called things. Following this direction, the Web-of-Things (WoT) paradigm has emerged with, as central vision, a common application layer targeting a seamless integration of heterogeneous devices [START_REF] Guinard | A Web of Things Application Architecture -Integrating the Real World into the Web[END_REF]. The strengths of this approach is to rely on well accepted protocols and patterns such as HTTP and Representational State Transfer (REST) services, providing high scalability while decoupling applications from software and hardware concerns.

We can nowadays find many of such smart things, for example sensors in buildings or daily-life objects in our homes. Smart environment applications are also emerging, aiming at increasing the user experience (or more generally the comfort) while optimising the energy consumption. Having access to historical data of sensors is often required by those applications. The field of smart buildings is a compelling example where control loops, or data-driven techniques targeting an intelligent adaptation of the building according to the users behaviour make an intensive use of historical data [START_REF] Ridi | Towards reliable stochastic data-driven models applied to the energy saving in buildings[END_REF]. Current approaches are deporting the storage outside of the sensor network in a single database or using a commercial cloud service. Although being viable solutions, they involve a series of constraints which are often not acceptable for building owners. The first worry is about the privacy of sensitive data such as presence sensor or electric consumption data. Secondly, forwarding the data outside the sensor network implies a gateway that is susceptible to be overloaded by the traffic and that represents a single point of failure as well as a point of entry for attacks. Considering those problematics, it appears that the best solution for owners would be to disconnect the sensor network from the rest of the world, thus increasing security while ensuring privacy by storing historical data within the sensor network.

In this paper, we present a mechanism for distributing historical sensor data on things, forming an in-network cloud relying on Web technologies. Our system distinguishes itself from others by being fully autonomous and requiring no human in the loop. We put special attention on limiting the network traffic while ensuring maximum fail safety and scalability.

RELATED WORK

Accessing data sources over RESTful APIs has made its way into databases especially built for storing time series. InfluxDB [1] is an open-source distributed databases making use of Web technologies such as HTTP, JSON and JavaScript for writing and reading data of sensors. Although it states itself as distributed and horizontally scalable, the documentation section describing the clustering is currently empty. StormDB [START_REF]Stormdb[END_REF] offers a cloud-based out of the box storage for time series. It bases on the Postgres-XC project for dis-tributing data across multiple Postgres databases. A load balancer dispatches requests to coordinators who are responsible for distributing those requests to data nodes. They are also responsible during read queries for merging the data coming from multiple data nodes. These freely available solutions only run on common operating systems as Linux, Mac and Windows. In addition, their core engines require a certain amount of resources (CPU and memory) that are not available on smart things. Applications using StormDB were indeed developed successfully for smart city context but relying on a remote cloud based implementation [START_REF] Difallah | Scalable anomaly detection for smart city infrastructure networks[END_REF].

When considering a distributed in-network storage, the placement of the data nodes is determinant when trying to minimize the total energy for gathering data to the storage nodes and replying queries. Several algorithms optimizing the placement of the data nodes are proposed in [START_REF] Sheng | Data storage placement in sensor networks[END_REF]. All of them are using trees as model structure for representing the problem. Trees are categorized as fixed or dynamic, depending on the problem's nature. Those algorithms, especially for dynamic trees are relevant to our problematic as devices can appear and disappear in a building.

The idea of storing historical data directly on sensor devices is not a new topic. This was studied for example in [START_REF] Mathur | Ultra-low power data storage for sensor networks[END_REF] where constrained devices used in sensor networks are considered for local storage. They developed Capsule, a rich, flexible and portable object storage abstraction that offers stream, file, array, queue and index storage objects for data storage and retrieval.

AIM AND ARCHITECTURE

We are proposing here an architecture, processes and interfaces supporting storage capacity on the already available devices in the sensor network, like sensors and gateways. With the advances made in electronics in recent times, sensors embed a non-negligible amount of memory that is often underused. This is even more the case for gateways having storage capacities of several mega-or even gigabytes. All this memory can therefore be used for storing historical data on Storage Peers (SP) that will be managed by Storage Coordinators (SC). We also propose to include the following features:

Location based grouping SP are placed among the network where data is produced to avoid traffic propagation over the whole network. This also ensures that a failure in a specific part of the building will not affect the rest of it. We hereby base our grouping according to the logical room structure of the building.

Replication for prevention An active agent takes actions for ensuring that data are replicated at least once on the SPs. Every SC is also backed up by a second one.

Energy efficient retrieval

The SC includes an algorithm based on the data repartition table to select the appropriate SP for limiting exchange of values.

Decoupling with network-layer The solution is decoupled from the network topology and routing protocol, able to work with IPv4 addresses as well with IPv6.

On-demand storage Storage resources are dynamically announced for saving storage space.

Efficient application layer CoAP (Constrained Application Protocol) [START_REF] Shelby | Constrained application protocol (coap)[END_REF] is used as application-layer protocol for accessing RESTful APIs. CoAP is lighter than HTTP and contributes to the efficiency of the system. The observe option [START_REF] Hartke | Observing resources in coap[END_REF] is used for implementing notifications between sensors and the storage peers.

Multicast based Multicast is used for efficient communication with several peers at one time in combination with the CoAP protocol.

Architecture

We propose an architecture based on the structural rooms and floors composition of a building. The logical structure of a building forms a tree where agents and producers are disposed. We extend this perception of a building by introducing the concept of zones grouping contiguous parts of a building, as illustrated in Figure 1. Each zone is managed by two SCs, one master and one backup, having their knowledge limited to their own content and to the references to their parent and children zones. This allows to split a building into autonomous regions that are not affected by failures in other zones. In addition, it avoids that each coordinator knows the whole structure of the building, therefore reducing the amount of synchronization data. Zones are formed dynamically during the operation according to the producers and the peers repartition. The principle of zone management is detailed in section 4.

!" !" !" !" # !# !# # # !# !# # # !# !# # !# !# # !# # !# # # !""#$%&"'()*"' *&+,-!(.$%/""!01%/""!2 !""#$%/""!2)&"'()*"' *&+,-!(.$%3 !""#$%/""!0)&"'()*"' *&+,-!(.$%,+4+.5 !""#$%,+4+.5)/""!0)&"'()*"' *&+,-!(.$%3

Entities and Roles

Our solution is composed of several agents having different behaviours depending mainly on their memory capacity. The agent applications are loaded in the device during deployment. Powerful devices may run multiple agents.

Storage Peers

Storage peers (SPs) manage the access to historical data according to the available local memory. They announce their capacity at startup to the Storage Coordinator (SC) managing the zone they reside in, indicating they are ready to receive data from producers. The storage management type used on the peer is hidden by a RESTful API providing functions for creating new datasets, adding entries as well as reading historical data. This approach allows for decoupling the storage from database technologies and exposes a single common storage interface for all peers. SCs and producers have no clue about which technology is used for the internal storage, as this could be for example MySQL, SQLite, flat files or even only EEPROM storage with specific indexing on constrained devices, for example using Capsule [START_REF] Mathur | Ultra-low power data storage for sensor networks[END_REF]. SPs regularly announce their storage capacity to the SC agent responsible for their location in the building.

Storage Coordinators

SCs manage the composition of the zones as well as the repartition of data between the peers. It starts up by looking after another coordinator already managing the location it resides in. If there is already a SC, it will put itself in stand-by, waiting to get promoted as backup. If no answer is received, it will promote itself as master coordinator for its location. It then periodically communicates with its parent and children in order to take decisions about merging or splitting zones. A major task of the storage coordinator is to manage the repartition of data between SPs by ensuring that each tuple of data is duplicated at least once. Finally, it serves as entry point for retrieval requests and will use its internal repartition table to find out which peers are concerned by the query.

Producer

A producer can be a sensor, a daily life object or any kind of thing that produces data in the form of notifications. The observer pattern of CoAP is here used among with multicast as notification mechanism for transmitting new values to storage peers in a single packet.

Client

The client stands for humans or applications requiring storage of historical data of a particular resource. Before a resource gets stored, an announcement by a client has to be sent to the storage coordinator responsible for the zone where the resource resides. Clients will afterwards be able to send read requests to the storage coordinator by indicating a resource and an interval of time for filtering purpose.

Interfaces

We define three types of interfaces that are used in our solution. Each one is accessible over a RESTful API and targets specific functionalities.

General management interface: SCs implement the server side of this multicast interface used for discovery purpose. Each participant of the cloud will use this interface to send a GET request for seeking after a SC managing a zone. For example, a client looking for the SC managing resources in the kitchen will send the following request: coap://229.0.5.32/home/floor1/kitchen. Each storage coordinator within the building will check if this location resides within their authority zone. Only the one that is managing this location will respond with the multicast address of its zone interface.

Zone interface: A specific zone interface is attributed to every zone, which is a unique multicast address generated during the zone creation. The purpose of this interface is for synchronization needs between the master coordinator with its backup, as well for the peers for communicating their filling ratio. In addition, this interface is used by the produc-ers for notifying peers with new values. We avoided CoAP groupcom [START_REF] Rahman | Group communication for coap[END_REF] because of constrained devices being unable to have multiple multicast addresses due to IP stack limitations. Clients send GET requests to this interface for retrieving historical data, as illustrated in the following request: coap://233.56.175.90/home/floor1/kitchen/temp?from= 2014-02-13&to=2014-02-18. Finally inter-zone synchronization is also achieved by using this interface.

Server interface: Each SC and SP has a server interface bound with its unicast address. For SCs, this interface allows to manage their participation in the observation of resources by giving the observe token. This interface also offers the possibility to retrieve the history of a resource stored on the peer by specifying an interval. The promotion of backup SCs is also managed using this interface.

The decomposition into several interfaces is an interesting approach, especially combined with multicast. First, having one interface for each behavior and particularly for each zone concentrates the communication within the specific zone without disturbing other zones or agents. Using multicast allows impersonal communications between participants. This impersonality improves the dynamics and fault tolerance of the solution as participants do not know the exact identity of their partners. We can illustrate this principle by citing the zone interface. SPs regularly notify their coordinator about remaining storage space. As this is achieved with multicast, they do not address a specific agent. This means that the agent can change over time without any consequence for the peers. For example if the master SC fails, peers will continue talking with the backup that will promote itself as master.

ZONE COMPOSITION STRATEGIES

As previously mentioned, our cloud storage is split into several zones forming a tree. Each zone is authoritative for one or many contiguous locations of the building. The location grouping process can be seen under several perspectives, and many algorithms can be imagined for defining zones. However, as our aim is also to build an energy efficient solution, we opted for an approach limiting the traffic between all the agents. SCs face two types of decisions during the grouping process: should they split their own zone into two new ones, and should they regroup a child zone with their current zone. Decisions are taken according to notions of costs.

Splitting

The coordinator for a zone Z will first compute the cost C(Z1, Z2) for each possible decomposition into two new zones Z1 and Z2 according to Eq. (1). We limit this computation to zones having at least two SCs and two SPs, and having obviously enough space for storing the historical data of resources residing within the zone. Finally, the minimum cost value among the possible decomposition is compared to the individual cost Ci(Z) of the zone Z. The zone will be split if the evaluation of Eq. 2 is true, i.e. if the best decomposition (minimising the cost) brings a gain of efficiency as compared to the current situation.

C(Z1, Z2) = Ci(Z1) + Ci(Z2) + Ce(Z1, Z2) (1) Ci(Z) -min(C(Z1, Z2)) > ǫs (2)
Ci Total number of hops (locations) between producers and storage peers within the zone. We consider that forwarding notifications to distant storage nodes is penalizing.

Ce Inter-zone costs depending on the filling rate of the child zone. Children zones are regularly sending some status data to their parent zone, which generates traffic. This traffic essentially depends on the filling rate of the child zone, as more messages will be exchanged as storage space decreases.

ǫs A cost gain factor above which the splitting is performed.

Merging

A SC can decide to merge its own zone with a child one in two situations:

Space: When one of the zones becomes full and has no more storage capacity. Grouping the zones will allow to increase the storage capacity and to postpone storage saturation.

Efficiency: When the cost of a single zone would be less than the actual situation. However, a SC has no knowledge about the internal structure (locations, storage peers and producers) of a child zone. This results in the incapacity of computing the costs for the merged zone. We opted for a solution where we compute an estimation of the cost Cp(Z) for the merged zone depending on each zone's costs and internal location tree depth, as described in Eq. ( 3). The merging is performed if the evaluation of Eq. ( 4) is true, i.e. if the merging brings a gain.

Cp(Z) = (Ci(Z1) + Ci(Z2)) * ln(D(Z1) + D(Z2)) (3)
C(Z1, Z2) -Cp(Z) > ǫm (4)
C(Z1, Z2) Two-zone cost computed as above.

D Depth of the zone's internal location tree. Each zone is composed of one or more contiguous locations represented as a tree. The depth of the tree influences the probability of having distant producers and consumers.

ǫm A cost gain factor above which the merging is performed.

PROTOCOL DETAILS

In this section we detail some behaviours of our proposed cloud storage for things. We assume that zones are already built and that the client has performed the SC discovery.

Data Announcement

We here describe the process of announcement and the internal communications as illustrated in Figure 2 part (a). The first step for a client is to announce its intention of storing data for a particular resource and for a given amount of time. This is realized by sending a multicast packet containing the necessary information to the SC. In order to offer a lightweight and robust mechanism for the announcement, we decided to specify the interchange format with a JSON schema, visible in Listing 1. JSON schemas [START_REF]Json schema[END_REF] are the equivalent to XML schemas and are used for validating data interchanged with Web services. As we are working in a world of constrained devices, using JSON appeared a reasonable choice.

Listing 1: JSON Schema for data storage announcement { "$schema ": "h t t p : / / j s o n -schema . o r g / d r a f t -04/schema #", " t i t l e ": " S t o r a g e " , " d e s c r i p t i o n ": "A r e s o u r c e s t o r a g e announcement " , "type ": " o b j e c t " , " p r o p e r t i e s ": { " u r l ": { " d e s c r i p t i o n ": "The r e s o u r c e t o s t o r e " , "type ": " u r i "} , "max-age ": { " d e s c r i p t i o n ": "How l o n g t h e data s h o u l d be s t o r e d " , "type ": " i n t e g e r "} , " u n i t ": { " d e s c r i p t i o n ": "Unit f o r max age " , "enum ": [ " day " , "month " , "y e a r "] } } , " r e q u i r e d ": [ " u r l " , "max-age " , " u n i t "] }

The SC performs a validation for ensuring that the received data is correct. It then ensures that the target resource is within its own zone. The selected appropriate SPs will be notified with a token that is used by CoAP for observing the resource. Finally, the SC enables the observation of the resource by providing the same token that was sent to the peers. In order to limit the network traffic, we combine observation and group communication. The resource will send the notification to the zone's multicast address, and only peers enabled for the token contained within the CoAP response will store the new data.

Data Retrieval

Once the announcement performed, clients have the possibility to retrieve the stored data, as depicted in Figure 2 part (b). This is achieved by the client sending a GET request to the SC responsible for the resource. The request must hold three parameters. The first one contains the URL of the resource, while the two remaining indicate the interval's start and end dates. This allows the client to filter the history by dates. The SC verifies that the requested resource is within its zone. If it is the case, it will then compute on which peers the data has to be retrieved in order to minimize the communications. Indeed, as each tuple of data is stored at least twice with overlapping intervals, this could result in unnecessary exchanges. The results are merged together in order to form a sorted JSON array of historical data. Due to the limited computational power of certain constrained devices, applying compression algorithms such as GZIP is not possible when exchanging historical data between SPs and SCs. However, JSON could be transposed in binary format for optimizing the communications.

EVALUATION AND DISCUSSION

In order to evaluate the performance of our system, and especially the retrieval of data, we have set an experimental zone. Our test zone was composed of one SC (Raspberry Pi) and three SPs (one Raspberry Pi and two OpenPicus Flyport). Each SP was preloaded with storage data as follows: 15072 entries for the Raspberry Pi, 14499 entries for the Flyport Wi-Fi and 14072 entries for the Flyport Ethernet. The SC was for its part preloaded with the repartition of data on the peers. To know the consequence of using a SC instead of directly talking to a REST data source, we set a second test bed only composed of a SP (Raspberry Pi). This peer was preloaded with the entire data set that was distributed among the peers in the previous configuration (without duplicate entries). The same retrieval request was sent for each test. It queries for historical data of a specific resource covering an interval of 12 hours. The answer for this request returns 39 entries, each one containing the timestamp of the measure and the associated value. Regarding the software stack, JCoAP 1 was updated to comply with the draft 18 of CoAP and used on the Raspberry Pi. For the Flyport, we opted for the currently lightest C implementation of CoAP, namely microcoap 2 . Both implementations of CoAP were improved for supporting multicast.

Distributed vs. Centralized

We compare both access types described above and show the round-trip time for each requests in Figure 3 part (a). First, each request is sent to the SC that will collect the distributed data on the SPs. For this case, the average round-trip time for 2500 consecutive requests is 186 milliseconds. In the second case, the client bypasses the SC and performs the requests directly on a storage peer storing the whole dataset. In this case, the average round-trip time was 70 milliseconds. Not surprisingly the performance for direct communication is better. We can deduce that the time used by the coordinator for retrieving the distributed data and building the response is about 115 milliseconds. However this additional delay is not that substantial so that the distributed approach would be limiting certain applica-1 https://code.google.com/p/jcoap/ 2 https://github.com/1248/microcoap tion scenarios.

Evaluating Concurrency

In this evaluation, we want to assess the difference between the approaches in terms of concurrent requests. Similarly to the previous case, we performed requests using the distributed approach and directly to the SP. Figure 3 part (b) shows the results when having up to 100 concurrent clients running 5 requests. Not astonishingly, the centralized approach scales much better as the rate of the distributed approach drops as soon as the number of concurrent clients is reaching 50. This limitation can be explained by the number of sub-requests that are sent by the SC to the peers. Indeed each client request will result in up to three peer requests. It is very likely that the SC, having limited memory and CPU is not able to handle such an amount of requests.

Discussion

Using our in-network cloud storage significantly eases the storage of historical data in a sensor network. No specific infrastructure is required for enabling the storage, as all available devices contribute in the storage effort by offering their capacities. While gateways offer large amount of memory, other things only have a few kilo-or megabytes at disposal. The total amount of memory within the network could thus be very limited and not enough for satisfying the requirements. In this case, we could imagine in the same way as for a cloud, to simply add more devices in the sensor network that will only be dedicated for storage, e.g. several Raspberry Pi offering cheap storage. The use of RESTful APIs on the peers storing time-series is an interesting alternative to the common SQL. Indeed smart things being mostly composed of constrained devices do not have the ability of implementing SQL. REST allows standardizing the way historical data are exposed in sensor networks, this in a very

!" !#$" !#%" !#&" !#'" !#(" !#)" !#*" %!!!"
%!(!" %$!!" %$(!" %%!!" %%(!" %&!!" %&(!" %'!!" %'(!" %(!!" +,-./01234"1356"789" +6:-68;"<" +,-./01234"1356"=646./3.>",.";?6"@;,2A>6"@;2A;6>B" =38;23C-;6/" D6.;2AE3F6/" (a)

!" #!" $!" %!" &!" '!" (!" )!" *!" +!" #!!" #" ##" $#" %#" &#" '#" (#" )#" *#" +#"
,-../00"123/"45"1/6-/030"789" :4;.-<</;3":=>/;30"?'"1/6-/030"@":=>/;3A" ,.2=2B>=>3C"D/03"E>3F":4;.-<</;3":=>/;30" G>03<>B-3/H" :/;3<2=>I/H" (b) simple manner, and probably sufficient for most scenarios.

The scalability test showed that the recomposition of data distributed on peers limits the number of concurrent clients. Nevertheless, it is unlikely that an average building management system issues such a high amount of history retrieval requests for a specific zone. Current machine-learning adaptation techniques only perform a few queries per hour, for which our proposal is more than compatible.

CONCLUSION

We have presented a system that composes a cloud-like storage space from several distributed smart things. Web technologies are at the heart of our system allowing a seamless integration into future building automation systems relying on the Web-of-Things as federating paradigm. The most relevant advantages of our approach are self-adaptation and autonomy. The storage automatically adapts itself to environmental changes when devices are added or removed by dynamically recomposing the zones. No human in the loop is required thanks to the agent technique that ensures auto-promotion and synchronization. Another important aspect is the efficiency that is part of every decision that the system takes. Traffic reduction is an important issue that is too often left aside by commercial technologies. Working in the field of smart buildings where the final aim is to save energy, it would be difficult to defend an additional layer that is itself not optimized.

Performance tests showed that even using very constrained devices as storage peers allows to obtain satisfying results in terms of round-trip time and scalability. Still, we plan optimizing the system in order to make it ready for future data-driven techniques requiring larger amount of historical data at higher frequencies. For this, we plan to compress the historical data exchanged between a storage coordinator and the peers, in order to accelerate the retrieval request. Furthermore, we intend to deploy our system for cloud storage in a real-life scenario to evaluate the splitting/merging strategy in a realistic situation.

Figure 1 :

 1 Figure 1: In-network cloud storage architecture based on building zones.

Figure 2 :

 2 Figure 2: The sequence of exchanged messages (a) for a client announcing data storage and (b) for a client retrieving historical values.

Figure 3 :

 3 Figure 3: Round-trip time for 500 consecutive requests (a) and scalability in terms of concurrent requests (b).

ACKNOWLEDGMENTS

The authors are grateful to the Swiss Hasler Foundation and to the RCSO grants from the HES-SO financing our research in this exciting area of smart buildings.