Régis Garnier

Christophe Guyeux
email: christophe.guyeux@univ-fcomte.fr

Stéphane Chrétien

Finding optimal finite biological sequences over finite alphabets: the OptiFin toolbox

In this paper, we present a toolbox for a specific optimization problem that frequently arises in bioinformatics or genomics. In this specific optimisation problem, the state space is a set of words of specified length over a finite alphabet. To each word is associated a score. The overall objective is to find the words which have the lowest possible score. This type of general optimization problem is encountered in e.g 3D conformation optimisation for protein structure prediction, or largest core genes subset discovery based on best supported phylogenetic tree for a set of species. In order to solve this problem, we propose a toolbox that can be easily launched using MPI and embeds 3 well-known metaheuristics. The toolbox is fully parametrized and well documented. It has been specifically designed to be easy modified and possibly improved by the user depending on the application, and does not require to be a computer scientist. We show that the toolbox performs very well on two difficult practical problems.

I. INTRODUCTION

Many biological data can be described as finite sequences over a finite alphabet. For instance, from a computer science viewpoint, a gene is nothing but a sequence of nucleotides A, C, G, T , and a protein is a sequence of amino acids (a.a., a set of 20 elements). In various situations, we can assign a score value to each finite sequence, that can account for its ability to solve a biological problem. For instance, these finite sequences can be the mapped to the descriptions of the protein's spatial configuration in the 3D square lattice; in this case, the score of the considered conformation is simply its energy. The main problem considered in this paper is the one of optimising the score value over all possible sequences. As is well known, most such problems are NP-hard and we often have to resort to efficient heuristics.

The main contribution of the present paper is to propose an handy and efficient toolbox that can solve such problems using a combination of metaheuristics. The toolbox is carefully designed in such a way that the sequences, their updates and the score function are easy to define and work with. Our toolbox -which is still in alpha stage, but in active development -incorporates the state of the art computation standards in computer science. It is mainly written in Python, it is MPI compliant, and it uses XML configuration files. The toolbox can be launched on a collection of computers through the network. After completion, the toolbox provides a report summarising various information concerning its convergence.

We provide examples of application to phylogenetics and protein structure predictions.

The remainder of this article is as follows. In the next section, we present the main problem and motivate its relevance to bioinformatics via some illustrative examples. In Section III, 3 currently embedded metaheuristics are presented. The toolbox's architecture is described in Section IV. Application examples are provided in Section V. The article ends with a conclusion, in which the contribution is summarized and future work is outlined.

II. WHAT TYPES OF PROBLEMS ?

A. Main objective Given an integer n > 0 and a finite set A, we consider the set S n of sequences of length n whose elements belong to A. In an equivalent manner, S n is the finite set of n-words w = w 1 ...w n whose letters w i belongs to the alphabet A. We consider a distance d : S n × S n -→ R + , and we suppose that d (S n , S n ) ⊂ N * : all distances between words are integer values. Finally, we consider a "scoring function" s : S n -→ R + , that is, a function that maps each n-word to a real positive value. Our objective is to design a software that uses various heuristics to find the minimum value of s on S n .

Another way of seeing is to construct an undirected graph. Each n-word can be identified with a node of a graph, while two nodes w and w are connected with an edge if and only if d(w, w ) = 1. At each node is associated its score. Our goal is then to design a way of moving from a node to the next in the graph, in such a way that the scoring function decreases over successive iterations. As will be clarified later, some sort of continuity will be needed : two nodes close enough must have, most of the times, approximately the same score. Some jumps will be acceptable, but should remain quite rare.

We will now motivate that this kind of problem arises in various instances in the field of bioinformatics by providing three important applications.

B. A first example: core genes and phylogeny

In the first example, we consider a set of species from which we want to infer a biomolecular phylogeny. For each species, we have a list of genes, i.e. a dictionary whose key is the gene name and value the DNA sequence. Although two given species do not necessarily have the same genes, any phylogenetic tree is based on the similarity of their common genes. Indeed, a phylogenetic tree can be computed as follows. Given a set of species, we extract a set of "core" genes (genes present everywhere in this set of species). Each core gene is multi-aligned using a toolbox such as e.g. Muscle [START_REF] Robert | Muscle: multiple sequence alignment with high accuracy and high throughput[END_REF], and then the alignments are concatenated. A mutation matrix between nucleotides is then built and the estimation of the tree corresponds to finding the most likely tree that, given this model of evolution, is able to explain the mutations found in the alignment. Statistical techniques like bootstrapping can provide support values on the obtained tree, and these supports provide information on the reliable character of the associated branch. In other words, a branch with a large support value will be considered as trustworthy.

For the same species, the phylogenetic trees inferred by using two different subsets of common genes are not equal in general. At least, it is unlikely that they share the same branch length. More dramatic, some sister relationship between species may be different (the trees do not have the same topology). Finally, the two obtained trees may share the same topology, but with different support values: some weakness may appear in the first tree in branches considered as trustworthy in the second one. If the given set of species have n genes in common (their core genome has n genes), then we can consider 2 n -1 different phylogenetic trees that show a certain variability in terms of topology and supports. To choose a tree that faithfully represents the phylogenetic relationship between the species at hand, we may consider that: (1) it is as well supported as possible, and (2) it is based on a nucleotidic sequence which is as large as possible. If such an approach is chosen, our objective is now to find the largest subset of core genes that leads to the tree that receives the largest support. This can be achieved as follows.

The n core genes are sorted alphabetically, and to each subset we associate a binary word of length n: its i-th character is 1 if and only if the i-th core gene is in the considered subset. The alphabet A is then equal to {0, 1}. Given a binary word w, after alignment, we can use a tool like RAxML [START_REF] Stamatakis | Raxmliii: a fast program for maximum likelihood-based inference of large phylogenetic trees[END_REF] to infer a tree as well as its supports. The score s(w) of this binary word can be the average value of the percentage of 1's inside w (number of considered core genes) and the lowest support within the tree (an integer ranging between 1 and 100). This provides us with a score we will have to maximise. We can also consider the inverse of this score, so as to turn the optimisation problem into a minimisation one.

C. A second example: self-avoiding walks and protein folding

Nowadays, it is relatively easy to have access to the genome of a given species, and to predict its coding sequences, either by blasting [START_REF] Altschul | Basic local alignment search tool[END_REF] a set of known genes (basis of knowledge) on it or by using ad hoc prediction tools like Glimmer for bacterias [START_REF] Arthur L Delcher | Identifying bacterial genes and endosymbiont dna with glimmer[END_REF]. It is thus quite easy to obtain the DNA sequence of a gene. A more difficult task is to predict its 3 dimensional shape, since this latter is mostly determined by the physico-chemical properties of its amino acid sequence, their coupling interactions, and so on. The protein prediction tools for conformation usually try to find the best conformation according to an optimization function. To do so, they focus on the hydrophylic/hydrophobic properties of each amino acid, and they try to find the 3D conformation that optimizes the number of hydrophobic a.a. inside the conformation (as proteins mainly evolve in aqueous solutions like cytoplasms).

In order to solve this optimization problem, a lattice is first considered, like the 2D or 3D square one, and then, a protein conformation in this lattice will be a self-avoiding walk (SAW [START_REF] Dill | Theory for the folding and stability of globular proteins[END_REF]), i.e. a succession of adjacent nodes of the lattice connected with edges, and such that each node strictly inside the walk has a degree of 2 (as two a.a. cannot lie in the same location), see Figures 1 and2. Given its sequence of hydrophobicity, finding the best 2D conformation of a protein is not an easy task. When considering the set of self-avoiding walks having n-steps and whose vertices are either black (hydrophobic) or white squares (hydrophylic residues), the authors of [START_REF] Crescenzi | On the complexity of protein folding (extended abstract)[END_REF] have indeed proven that determining the SAWs that maximize the number of neighboring black squares in this set, which is supposed to be the conformation chosen by the protein, is NP-hard.

Given a sequence of amino acids, we will resort to using relevant heuristics for predicting the most probable conformation of the protein. The translation in terms of S n and scoring function of Section II-A is obvious. Each self-avoiding walk is a word on the {North, East, South, West} alphabet, which describes its absolute encoding. The scoring function, here, is the number of neighboring black squares.

D. A third example regarding the protein folding process

The problem described in the previous subsection also leads naturally to the study of various subsets of self-avoiding walks [START_REF] Bahi | Protein folding in the 2D hydrophobic-hydrophilic (HP) square lattice model is chaotic[END_REF], [START_REF] Bahi | Chaos of protein folding[END_REF], [START_REF] Guyeux | Is protein folding problem really a NP-complete one? first investigations[END_REF]. In the first approach, starting from the straight line, we obtain by a succession of pivot moves of ±90 • (to choose a node as center of a rotation of the tail) a final conformation being a self-avoiding walk. All the intermediate conformations must be self-avoiding too (see Fig. 1). Such a procedure is one of the two most standard implementations of the so-called "SAW requirement" in the bioinformatics literature, which is for instance presented in [START_REF] Islam | Clustered memetic algorithm for protein structure prediction[END_REF], [START_REF] Unger | Genetic algorithm for 3d protein folding simulations[END_REF], [START_REF] Higgs | Genetic algorithm feature-based resampling for protein structure prediction[END_REF]. We have previously shown in [START_REF] Guyeux | Is protein folding problem really a NP-complete one? first investigations[END_REF] that, in this first category of prediction software, it is impossible to reach all the possible conformations.

The other main approach of protein structure prediction using self-avoiding walks starts with an 1-step SAW, and at iteration k, a new branch is appended to the current tail of the walk, in such a way that the new k-step self-avoiding walk has a larger number of neighboring black squares (see Fig 2). The protein is thus constructed step by step, reaching the best local conformation at each iteration. It is easy to see that such an approach leads to all the possible self-avoiding walks having the length of the considered protein [START_REF] Guyeux | Is protein folding problem really a NP-complete one? first investigations[END_REF].

As stated above, these two approaches cannot reach the same set of conformations. One very bad problem faced by the community is that the difference between the two sets is far from negligible, leading to the necessity of comparing their respective sizes. Clearly, non unfoldable self-avoiding walks, that is, SAWs on which any pivot move introduces an intersection between the head and the tail of the new walk (like in Fig. 3), cannot be reached as the result of performing iterative pivot moves starting with the straight line. Our goal is to generate unfoldable self-avoiding walks. In order to do this, a relevant score will be the number of possible pivot moves, which must reduce to zero in the case of an unfoldable selfavoiding walk.

III. THREE EMBEDDED HEURISTICS

In the present section, we outline the three heuristics that we have embedded within our toolbox. Further details can be found in [START_REF] Alkindy | Hybrid genetic algorithm and lasso test approach for inferring well supported phylogenetic trees based on subsets of chloroplastic core genes[END_REF], [START_REF] Alsrraj | Well-supported phylogenies using largest subsets of core-genes by discrete particle swarm optimization[END_REF] or in the software documentation 1 .

A. Genetic algorithm approach

The population is initialized with a collection of 50 random words. Then, the genetic algorithm iterates until discovering a word whose score is larger than a threshold, or at most for 200 iterations. Each iteration, which produces a new population, consists of the following steps:

1) Repeat 5 times a random pickup of a pair (w 1 , w 2 ) of words and mix them using a crossover approach. In this step, indexes {1, . . . , n} are partitioned into k, k ≤ n 2 , subsets I 1 , . . . I k . A new word w is then defined by w i = w 1 i if i belongs to some I j where j is odd; otherwise w i = w 2 i . The obtained words are added to the population P , resulting in population P c . 2) Mutate 5 words of the population P c . More precisely, for each of these words w, k randomly selected values of w are replaced by a character randomly picked in A, leading to a new word. The mutated words are added to P c leading to population P m . 3) Produce population P r by adding 5 new random words to P m . 4) Select the 50 best words in population P r to form the new population P .

The whole set of produced words with their associated scores contains valuable information about which parts of the words tend to decrease the scores. A Lasso test [START_REF] Tibshirani | Regression shrinkage and selection via the lasso[END_REF] may then be applied, to rank the word positions according to their impact on the scoring function. In this case, a last genetic algorithm phase is launched on the updated population, in order to mix these candidate words.

B. Particle swarm optimization approach

The particle position is a vector of N values belonging in A, an alphabet of positive numbers. The objective is to define a way to move the particles in the N dimensional search space so that they produce the optimal vectors with respect to the scoring function. A swarm of L particles is then a list of position vectors (X 1 , X 2 , . . . , X L ) together with their associated velocities (V 1 , V 2 , ..., V L ), which are n-dimensional vectors of real numbers between 0 and 1. The latter are initialized randomly. At each iteration, a new velocity vector is computed as follows:

V i (t + 1) = wV i (t) + φ 1 P best i -X i + φ 2 P best g -X i , (1)
where w, φ 1 , and φ 2 are weighted parameters setting the level of each three trends for the particle, which are respectively: to continue in its adventurous direction, to move in the direction of its own best position P best i , or to follow the gregarious instinct to the global best known solution P best g . Both P best i and P best g are computed according to the scoring function. The new position of the particle is then obtained, by applying the sigmoid function [START_REF] Mojtaba | A novel binary particle swarm optimization[END_REF] on each velocity, and by considering the value of A the most close to this latter.

C. Simulated annealing approach

After an initialization step, the Simulated Annealing loop is composed by (a) a move in the neighborhood of the current solution, (b) an evaluation of this new position by a real-valued scoring function, then (c) a test, given a well chosen criterion, to store this position as the new best one. Various criteria can be considered, they all have been embedded within the software. A detailed description of the simulated annealing algorithm is provided in Figure 4.

We are now able to describe the software that embeds these metaheuristics, in order to solve the general problem of Section II-A.

IV. SOFTWARE DESCRIPTION

A. Requirement Guidelines

Let us first outline the requirements that were set up before the conception of our toolbox.

• We need to assess the efficiency of the meta-heuristics (MH) for various models of solution spaces. Therefore, we focus on measuring and keep track of the performances, record events or messages (logfiles) related to our algorithms. In other words, introspection is preferred to speed. • We need to allow for easy application to new problems the users might come up with. Similarly, new algorithms and heuristics must be easy to embed, as well as batteries of tests. • As we intend to work on multi-cores or clusters and with scoring functions that putatively need a long computation time, we look for a scalable, ready, and reliable network solution. Furthermore, as the heuristic methods under consideration need various amount of resources during computation, we need a convenient scalability mechanism. For this purpose, we chose the Message Passing Interface (MPI) 2.2. • Various existing scripts and third party applications need to be embedded or reused without code duplication, and both must be launchable for testing purposes or specific needs. Therefore, we need an easy to use but reliable shell script launcher, with basic communication and logging capabilities. • Finally, targeted end-users like bioinformaticians are typically familiar with the use of terminals and code scripting, but they are not necessarily computer scientists. This is why the Python language has been chosen. But with a Clike flavor syntax when possible, to produce codes easy to read for everyone.

T ← calc(T) s0 (

B. Architecture Guideline

Software architecture is meant to convert requirements in a small list of well-known design patterns. Objectives of the latter is to provide a coherent approach that allows to easily reuse, maintain, and upgrade the proposal. Indeed, a pattern is used whenever possible, while they are reused as soon as possible to reduce their number, leading to an architecture of low complexity as depicted in Fig. 5. It is described hereafter. 1) Introspection -Decorators: Oriented Object Programming (OOP) is the main chosen paradigm of the framework but, as Python offers convenient ways to use it, a kind of Aspect Oriented Programming (AOP) has been preferred to separate some technical functions from domain ones. Python decorator (wrapper) has been the chosen pattern, which is in this language a simple decoration on a class, function, or method. Basically a decorator provides a way to intercept communication to or from a designated function. So inputs, outputs, and capability of a decorated class could be analyzed (loggers), validated (contracts), have methods replaced (monkey patching) or enhanced without modifying its code. Without the necessity to launch tests again on the decorated domain specific code. The function domain is also more readable and simpler to write. Note that a decorator can be shared with as much domain functions as needed. So decorators maximise the reuse of codes, leading to a decrease of both the testing time and the failure risk.

2) Framework -low coupling, helpers, and DAO: a) Coupling and cohesion, Liskov principle: Coupling property is a quality metric that describes the interdependence degree between software modules. Cohesion, for its part, is another metric that describes how much elements of a given module are related on. Our objective was to increase as much as possible these two metrics, for each module within the software. As all MJs use both move and score concepts (we called them "field functions"), separated libraries have been built for these functions. Indeed, the generalization of a concrete move function leads to a common interface with all other move functions. Some characteristics are also shared by all score functions, which leads to a second interface. So, each member of our class of functions respects the so-called Liskov substitution principle: if S is a subtype of T, then objects of type T can be replaced with objects of type S. b) Factory, generic programming, and mapping: If we have substitutive functions but call a specific type of them in the main MH template, we need to change code of template each time we change experiment conditions. The coupling must be lower. We could call a function that is interface compatible. So we need a tool to select function in a cohesive collection from a XML hashtables and return instance of this function. This is provided by the factory and dependency injection pattern.

As no reference to effective calculus function remains in the main MH code, we abusively speak about templates, in a C++ terminology. Launcher scripts are mostly strategy patterns, i.e., they could be launched by any elemnts of the templates. By doing so, the targeted quality metric was genericity: it measures the ability for effective types of data to be defined dynamically, i.e., only when instantiated. With this low-decoupling approach, the framework allows a zeromodification reuse of old field functions with a new template, and vice versa. c) Template helpers: MH templates often share same concepts but with specific implementation. Associated functions of common concept are outsourced. This approach has been chosen as we wanted not to allow any lambda user to practice a top-down programming. Furthermore, the objective was to discourage a "scripting" bottom up style for templates, as authors considered that users will want to plan neither unit nor integration tests on templates. By doing so, code review is more easy to achieve. d) DAO and loggers: As stated previously, various formerly existing shell scripts were at the beginning of this software, and we tried to group them in a unique benchmarkable and reliable environment. Scripts usually pipe applications, while inputs and outputs of these applications are usually files. To be readable, all these files must follow specific standards. As a consequence, conversions, sorting, or subset selection are sometimes required regarding these files. To prevent from issues that may be raised by such a variety of inputs and outputs, two priorities have been kept in mind: an improved and stringent logging and a file system (Data Access Object, i.e., DAO pattern) helpers.

3) Maintainability -Third-party libraries policy: Authors of this software have a priori considered that Python is not the most used language, and so developers may not be familiar with its ecosystem. This lack of background does not allow all developers to sort available libraries according to their maturity or their reliability, while various libraries from the Python ecosystem are still in an early reversion designed from scratch. Consequently, we decided to use as few dependencies as possible, leading to a very simple proposal easy to improve, customize, and fork.

4) Scalability -MPI 2.2 and sockets libraries, server MPI factory: Let us firstly recall that Message Passing Interface (MPI) is a concurrency inter-process communication system based on message-passing. The message-passing paradigm implies than process code execution is message-driven, i.e., the type and content of message select the code to invoke for its receiver. MPI is a Language Independent Specifications (LIS), point to point, and multicast communication protocol. It mainly provides a set of basic tools for synchronisation and data exchange, being an ISO 5th, 6th, and 7th level middleware and application framework using TCP sockets. APIs are provided for C, C++, and Fortran languages so, in Python, a binding library or a MPI standard implementation is needed. This latter is described hereafter.

a) MPI4py and OpenMPI: MPI for Python (MPI4py) is based on the standard MPI 2.2 C++ bindings and it provides a similar interface. But MPI4py adds an unique Python capability, pickling, to send/receive objects in RAM or HDD buffers. The widely used middleware OpenMPI 2 has been chosen, as this is a software of production quality which is fully MPI 2.2 and 3.1 compliant. It allows process spawning, network and process fault tolerance, safe threads, a good portability, and so on. Last, but not least, it is well documented and supported.

b) Client/server mode and libraries: MPI4py over Open-MPI third party library was a prerequisite as most of existing scripts use it. However, after an intensive evaluation, its low reliably with long series of connections, spawning, or big message exchanges was quite disappointing. Unfortunately, no good alternative has been found.

As a consequence, we decided to focus on the MPI 2.2 client/server (CS) mode, while minimizing as possible the coupling between MPI4py and the other codes within the software. Indeed, only a MPI-helper library requires it. c) MPI process server factory: Using OpenMPI name server, a socket helper library and the OpenMPI helper, we allow a MPI process to launch another MPI 2.2 subprocess (Fig. 7). And to free resources after results return. This server XML configured factory with socket communication channel is a standalone script (singleton).

d) MPI most used pattern -Master and Workers: In this parallel computation pattern, a main thread called the master splits a given problem and merges partial results effectively computed by a set of forked processes, which are called the workers. This pattern can rely on a messagepassing architecture to broadcast data to workers and regroup partial obtained results. This is typically a good choice when a problem can be sliced in a set of same smaller and simpler ones, splitting and merging times being small when compared with sub-problems computation ones.

5) Reuse with low-coupling -Executor library: The MPI server factory needs to execute MPI and Python codes in a shell, and with a maximum of security. As a set of applications and scripts has been added that use XML configuration files for their environment, a shell launcher has been required, which should follow an "easy to use" philosophy like in Java 5 executor with easy contracting and centralized logging abilities.

6) Introspection -general architecture and dataflow programming: According to the observation granularity, the proposed framework uses various classical software architectures. Some components as the MPI factory server induce their own main architecture, and the classical architecture erosion (gap between models and implementation) is observed. This is why, as much as possible, Separation of Concerns (SoC) is used at high level structures.

But a significant choice has been forced and preserved over versions concerning the "state" data of each MH process. Each problem (that is, the combination of a MH script, field functions like the move or score ones, and their parameters) is defined in an unique XML file. This XML is validated by XSD; it is parsed, and then converted in a Python dictionary, i.e., a hash-table implementation. So (almost) all the needs start in an unique entity. The choice has been made to preserve this entity, and to add in it all what MH templates require to be instantiated. In other words, this is principally a memory for current coordinates and score when moving on the solution space. This state token, which is a dataflow programminginspired architecture, is passed to each function and all of them return new version of the token.

With this approach, MH templates can be modeled by graphs, thus simplifying both software monitoring and tuning. Additionally, the check-pointing strategy that corresponds to store this token leads to various interesting properties, like the possibility to restore a computation from a given saved point, to study the convergence speed, to benchmark various physical platforms on a same computation, or to launch in parallel various instances to maximize the chance to find interesting solutions (or to study the dispersion of solutions). We thus improve both introspection and reliability.

C. How to construct a MH in practice

As a framework, a toolbox for computing resources is provided, together with a methodology to create MH. In broad outline some simple templates are created, helper are written with DOA and helpers, then when needed, MPI processes are nested with MPI helper and standalone MPI C/S factory. Figure 8 illustrates a SA evolution build by iterating SA and greedy local descent basic algorithms. The first step was to construct both SA and local descent in an one-core version. They are then linked by a constructor to apply a descent on each solution selected by the SA Metropolis criterion. In a second step, the descent is replaced by a MPI version: some descents are launched in Workers and the Master returns the best Worker result. The simple Executor is replaced by a call to the MPI server factory. Then, in a third step, the simple SA execution itself is injected in a Worker and a Master, created to make a synchronisation after each Markov chain.

On Figure 8, dashed arrows represent data circulation, and then moves of the token instances.

V. APPLICATIONS

A. Protein folding

The first illustrative example of the effectiveness of the proposed software is related to the protein folding problem presented in Section II-D. More specifically, given a sequence length, the objective is to find a self-avoiding walk of this length which is non unfoldable: to apply a pivot move at each internal node always leads to a new walk that does not satisfy the self-avoiding requirement. We considered the number of nodes on which a pivot move can be applied as the scoring function. When this score is zero, then a non unfoldable self-avoiding walk has been obtained. As a first basic experiment, we considered the SA metaheuristic on short sequences, to reduce the computation time. Obtained results are depicted in Figure 9, which emphasize the convergence ability of the proposed software on the considered problem.

B. Phylogeny

We have first considered phylogenetic issues, with the family listed in Table I. This family is constituted by 6 genomes, with a number of detected genes that ranges from 138 to 271, and a core genome of 122. The phylogeny with the alignment of these core genes leads to a small weakness in one branch (bootstrap lower than 95). To wonder whether some genes may be responsible of such weak uncertainty, we have firstly launched the genetic algorithm, which has stopped after 29 iterations, in systematic investigation mode, leading to 2 topologies:

• Topology 0 depicted in Fig. V-B, that has occurred 27 times. The best obtained tree has a lowest bootstrap of 96, while in average the lowest bootstrap is equal to 86. The PSO, for its part, has rapidly found a first well supported phylogenetic tree in a third different topology, and with all supports equal to 100, namely: (((C.closterium,F.solaris),L.danicus),C.daemon). However, the PSO has used only 47.5% of the core genes to reach such a tree. According to our stop criterion, this tree has not been returned by the algorithm. Indeed, this example illustrates the ability of the particle swarm optimization algorithm to more globally visit the whole space at the beginning, in order to discover regions of interest.

The simulated annealing, for its part, raised these 3 topologies. A remarkable element is that these 3 topologies have the whole bootstraps equal to 100. Furthermore, Topology 2 appears as the best one according to the produced result (it was Topology 0 according to the GA, while PSO has not succeeded in separating these two topologies). Obviously, both PSO and SA have converged to local minima that are not global ones if we consider that both minimum bootstraps and proportion of core genes must be maximized. Launching them again with other initial values and parameters may select other optimal positions in the cube.

VI. CONCLUSION

In this article, we have detailed the general principles at the bases of an efficient software we propose for the bioinformatician community. This software iterates on finite words on a finite alphabet, such that each word has a score. Using one of the three most classical metaheuristics, this software is able to discover the optimums in a distributed manner. The software is easy to configure, to use, and to update. It can resolve a variety of bioinformaticians problems, like to measure the effects of genes on a phylogeny.

In future work, we intend to improve the following elements. First of all, the XML "problem" configuration file must have a more common structure, to allow state token to be used by more than one specific MH template. The objective is as follows: a problem, started to be resolved by a SA, should be finalized by a PSO, and vice versa. An improvement of the execution speed must be achieved, among other things by substituting some of our functions by more efficient and reliable ones provided by third-party libraries. Finally, the number of tests, contracts, and code factoring must be enlarged, together with logs, error management, and documentation, to help bug tracking, maintainability, safety, and data integrity.

Fig. 1 .

 1 Fig. 1. Protein Structure Prediction by folding SAWs

Fig. 2 .

 2 Fig. 2. Protein Structure Prediction by stretching SAWs

Fig. 3 .

 3 Fig. 3. An unfoldable self-avoiding walk

Fig. 4 .

 4 Fig. 4. Simulated annealing as a threshold class algorithm.

Fig. 5 .

 5 Fig. 5. Detailed architecture

Fig. 6 .

 6 Fig. 6. Dependency Injection, Factory and Mapping Example (slightly simplified)

Fig. 7 .

 7 Fig. 7. Multi-level MPI 2.2 with Server MPI Factory and Executor Class functioning Principles (slightly simplified)

Fig. 8 .

 8 Fig. 8. 2-level MPI 2.2 application : SA with Basin Hopping algorithm (simplified)

Fig. 9 .

 9 Fig. 9. Example of output exploitation using the SA on a short size problem.

Fig. 10 .

 10 Fig. 10. First obtained supported topology.

The software is currently downloadable on demand at https://bitbucket. org/rjpGarnier/projetm2recuit.git


Abdessalem Abidi
email: abdessalemabidi9@gmail.com

Christophe Guyeux

Belgacem Bouallègue

Mohsen Machhout

Conditions to have a well-disordered dynamics in the CBC Mode of Operation

published or not. The documents may come

I. INTRODUCTION

Block ciphers, like the Data Encryption Standard (DES) or the Advanced Encryption Standard (AES), have a very simple principle. They do not treat the original text bit by bit but they manipulate blocks of text. So, it is not sufficient to put anyhow a block cipher algorithm in a program. We can, instead, use these algorithms in various ways according to their specific needs. These ways are called the block cipher modes of operation. There are several modes of operation and each one possesses its own characteristics and its specific security properties. In this article, we are only interested in one of these modes, namely the cipher block chaining (CBC) mode, and we will study its dynamical behavior of chaos.

The chaos theory that has been consider in this article is the Devaney's topological one and its ramifications [START_REF] Devaney | An Introduction to Chaotic Dynamical Systems[END_REF]. Being reputed as one of the best mathematical definition of chaos, this theory offers a framework with qualitative and quantitative tools to evaluate the notion of unpre-dictability [START_REF] Jacques M Bahi | Efficient and cryptographically secure generation of chaotic pseudorandom numbers on gpu[END_REF]. As an application of our fundamental results, It has been interested in the area of information safety and security. Specifically, the contribution belongs to the field of the cipher block chaining modes of operation.

In [START_REF] Abidi | Proving chaotic behavior of cbc mode of operation[END_REF], we have been started to give mathematical proofs that emphasize the chaotic behavior of the CBC mode of operation. Thereafter, in [START_REF] Abidi | Quantitative evaluation of chaotic cbc mode of operation[END_REF], we have been stated that in addition to being chaotic as defined in Devaney's formulation, this mode is indeed largely sensible to initial errors or modification on either the Initialization Vector IV or the message to encrypt. Its expansivity has been regarded too, but this property is not satisfied as it has been established thanks to a counter example. In this new article, we have been intend to deepen the topological study of this CBC mode of operation in order to obtain a complete mathematical overview of its dynamics.

The remainder of this research work is organized as follows. In the next section, some basic definitions related to chaos and cipher block chaining mode of operation will be recalled. Previously obtained results are recalled in Section III. Sections IV and V contain the main contribution of this article: the first one evaluates the topological mixing of the CBC mode of operation, while the second one focuses on its topological entropy. This article ends with a conclusion section where the contribution is summarized and intended future work is presented.

II. BASIC RECALLS

This section is devoted to basic definitions and terminologies in the field of topological chaos and in the one of block cipher mode of operation.

A. Devaney's Chaotic Dynamical Systems

In the remainder of this article, m n denotes the n th block message of a sequence S while m j stands for the j -th bit of integer of the block message m ∈ 0, 2 N -1 , expressed in the binary numeral system and x i stands for the i th component of a vector x.

X N is the set of all sequences whose elements belong to X .

f •k = f •...•f is for the k th composition of a function f . N is the set of natural (non-negative) numbers, while N * stands for the positive integers 1, 2, 3, . . . Finally, the following notation is used:

1; N = {1, 2, . . . , N }.
Consider a topological space (X , τ ), where τ represents a family of subsets of X , and a continuous function f : X → X on (X , τ ).

Definition 1

The function f is topologically transitive if, for any pair of nonempty open sets U, V ⊂ X , there exists an integer k > 0 such that f

•k (U) ∩ V = ∅. Definition 2 An element x is a periodic point for f of period n ∈ N, n > 1, if f •n (x) = x and f •k (x) = x, 1 ≤ k ≤ n.
Definition 3 f is regular on (X , τ ) if the set of periodic points for f is dense in X : for any point x in X , any neighborhood of x contains at least one periodic point.

Definition 4

The function f has sensitive dependence on initial conditions on the metric space (X , d) if there exists δ > 0 such that, for any x ∈ X and any neighborhood V of x, there exist y ∈ V and n > 0 such that the distance d between the results of their n th composition, f •n (x) and f •n (y), is greater than δ:

d (f •n (x), f •n (y)) > δ.
δ is called the constant of sensitivity of f .

Definition 5 (Devaney's formulation of chaos [1])

The function f is chaotic on a metric space (X , d) if f is regular, topologically transitive, and has sensitive dependence on initial conditions. Banks et al. have proven in [START_REF] Banks | On devaney's definition of chaos[END_REF] that when f is regular and transitive on a metric space (X , d), then f has the property of sensitive dependence on initial conditions. This is why chaos can be formulated too in a topological space (X , τ ): in that situation, chaos is obtained when f is regular and topologically transitive. Note that the transitivity property is often obtained as a consequence of the strong transitivity one, which is defined below. Definition 6 f is strongly transitive on (X , d) if, for all point x, y ∈ X and for all neighborhood V of x, it exists n ∈ N and x ∈ V such that f •n (x ) = y.

B. CBC properties

As what has been already defined, a mode of operation is an algorithm that uses a block cipher to provide an information service such as confidentiality or authenticity. The most commonly used mode of operation in recent decades is the cipher block chaining CBC. In what follows, it will be seen how this mode works in practice.

1) Initialization vector IV: Like some other modes of operation, the CBC one requires not only a plaintext but also an initialization vector (denoted as IV in what follows). An IV is an arbitrary number that must be generated for each execution of the encryption operation. For the decryption algorithm, the vector used must be the same, see Figure 1. The use of this vector prevents repetition in data encryption. So, it offers the benefit to make this operation more difficult for a hacker. Indeed, it deprives him to find patterns and break a cipher. The length of the initialization vector (the number of bits or bytes it contains) depends on the method of encryption. It is usually comparable to the length of the encryption key or block of the cipher in use.

2) CBC mode characteristics: In cryptography, Cipher Block Chaining is a block cipher mode that provides confidentiality but not message integrity. It offers a solution to most of the problems presented by the Electronic Code Book (ECB) mode; thanks to CBC mode, the encryption will depend not only on the plaintext, but also on all preceding blocks. More precisely, each block of plaintext is XORed immediately with the previous cipher text block before being encrypted (i.e., the binary operator XOR is applied between two stated blocks). For the first block, the initialization vector acts as the previous cipher text block, see Figure 1. CBC mode possesses several advantages. In fact, the same plaintext is encrypted differently in the case of different initialization vectors. In addition, the encryption operation of each block depends on the preceding one, so any modification in the order of the cipher text block makes the decryption operation unrealizable. Furthermore, if a transmission error affects the encrypted data, for example the C i block, so only m i and m i+1 blocks will be influenced by this error while the remained blocks will be determined correctly. Conversely, the CBC mode is characterized by two main drawbacks. The first one is that encryption is sequential (i.e., it cannot be parallelized). The second one is that the initial message should be padded to a multiple of the cipher block size.

In the next section, the previous results that has been detailed respectively in [START_REF] Abidi | Proving chaotic behavior of cbc mode of operation[END_REF] and [START_REF] Abidi | Quantitative evaluation of chaotic cbc mode of operation[END_REF] will be summarized.

A. Modeling the CBC mode as a dynamical system

The modeling follows a same canvas as what has been done for hash functions [START_REF] Bahi | Hash functions using chaotic iterations[END_REF], [START_REF] Guyeux | A topological study of chaotic iterations. application to hash functions[END_REF] or pseudorandom number generation [START_REF] Bahi | Evaluating quality of chaotic pseudo-random generators. application to information hiding[END_REF].

Let us consider the CBC mode of operation with a keyed encryption function E κ : B N → B N depending on a secret key κ, where N is the size for the block cipher, and D κ : B N → B N is the associated decryption function, which is such that ∀κ, E κ • D κ is the identity function. We define the Cartesian product X = B N ×S N , where:

• B = {0, 1} is the set of Boolean values, • S N = 0, 2 N -1 N ,
the set of infinite sequences of natural integers bounded by 2 N -1, or the set of infinite N-bits block messages, in such a way that X is constituted by couples: the internal states of the mode of operation, and sequences of block messages. Let us consider the initial function:

ι : S N -→ 0, 2 N -1 (m i ) i∈N -→ m 0
that returns the first block of a (infinite) message, and the shift function:

σ : S N -→ S N (m 0 , m 1 , m 2 , ...) -→ (m 1 , m 2 , m 3 , ...)
that removes the first block of a message, when counting from the left. We define:

F f : B N × 0, 2 N -1 -→ B N (x, m) -→ x j m j + f (x) j m j j=1..N .
This function returns the inputted binary vector x, whose m j -th components x m j have been replaced by f (x) m j , for all j = 1..N such that m j = 0. In case where f is the vectorial negation, this function will correspond to one XOR between the plaintext and the previous encrypted state. The CBC mode of operation can be rewritten as the following dynamical system:

X 0 = (IV, m) X n+1 = (E κ • F f0 (ι(X n 1 ), X n 2 ) , σ(X n 1 )) (1)
For any given g : 0, 2 N -1 × B N -→ B N , we denote G g (X) = (g(ι(X 1 ), X 2 ); σ(X 1 )) (when g = E κ • F f0 , we obtain one cipher block of the CBC, as depicted in Figure 1). The recurrent relation of Eq.( 1) can be rewritten in a condensed way, as follows.

X n+1 = G Eκ•F f 0 (X n ) . (2)
With such a rewriting, one iterate of the discrete dynamical system above corresponds exactly to one cipher block in the CBC mode of operation. Note that the second component of this system is a subshift of finite type that is related to the symbolic dynamical systems known for their relation with chaos [START_REF] Lind | An introduction to symbolic dynamics and coding[END_REF]. Now, a distance on X as follows: d((x, m); (x, m)) = d e (x, x) + d m (m, m) is defined, where:

         d e (x, x) = N k=1 |x k -xk | d m (m, m) = 9 N ∞ k=1 N i=1 m i -mi 10 k .
This distance has been introduced to satisfy the following requirements:

binary components that are different between the two internal states X 1 and Y 1 . • The k-th digit in the decimal part of the distance between X and Y is equal to 0 if and only if the k-th blocks of messages X 2 and Y 2 are equal. This desire is at the origin of the normalization factor 9 N .

B. Proofs of chaos

As mentioned in Definition 5, a function f is chaotic on (X , τ ) if f is regular and topologically transitive. It has been beagan in [10] by stating some propositions that are primarily required in order to proof the chaotic behavior of the CBC mode of operation.

Proposition 1 Let g = E κ • F f0
, where E κ is a given keyed block cipher and f 0 : B N -→ B N , (x 1 , ..., x N ) -→ (x 1 , ..., x N ) is the Boolean vectorial negation. The directed graph G g is considered, where:

• vertices are all the N-bit words.

• there is an edge m ∈ 0, 2 N -1 from x to x if and only if g(m, x) = x. If G g is strongly connected, then G g is strongly transitive.

Then it has been proven that, Proposition 2 If G g is strongly connected, then G g is regular.

According to Propositions 1 and 2, It can be conclude that, depending on g, if the directed graph G g is strongly connected, then the CBC mode of operation is chaotic according to Devaney, as established in our previous research work [START_REF] Abidi | Proving chaotic behavior of cbc mode of operation[END_REF]. In this article and for illustration purpose, we have been also given some examples of encryption functions making this mode a chaotic one.

In the next section, some quantitative measures of chaos that have already been proven in previous research work will be recalled.

C. Quantitatives measures

In [START_REF] Abidi | Quantitative evaluation of chaotic cbc mode of operation[END_REF], these two following propositions have been respectively developed .

Proposition 3

The CBC mode of operation is sensible to the initial condition, and its constant of sensibility is larger than the length N of the block size.

Proposition 4

The CBC mode of operation is not expansive.

To sum up, CBC mode of operation is sensible to the initial conditions but it is not expansive.

IV. TOPOLOGICAL MIXING

The topological mixing is a strong version of transitivity.

Definition 7 A discrete dynamical system is said topologically mixing if and only if, for any couple of disjoint open set U, V = ∅, there exists an integer n 0 ∈ N such that, for all n > n 0 , f

•n (U) ∩ V = ∅. Proposition 5 (X , G g ) is topologically mixing.
This result is an immediate consequence of the lemma below.

Lemma 1

For any open ball B = B((x, m), ε) of X , an index n can be found such that G •n g (B) = X . PROOF Let B be an open ball whose radius ε can be considered as strictly lower than 1, ε < 1.

All the elements of B have the same state x, and due to the definition of the chosen metric, they are such that an integer k (= -log 10 (ε) ) satisfies:

• all the strategies of B have the same k first block messages, • after the index k, all values are possible. Then, after k iterations, the new state of the system is G •k g (x, m) 1 and all the strategies are possibles (any point of the form

(G •k g (x, m) 1 , m), with any m ∈ S, is reachable from B).
Let (x , m ) ∈ X . We will prove that it can be reached by starting from B. Indeed, let us consider the point (x, m) of B defined by:

• x = x • ∀i k, mi = m i , • mk+1 = G •k g ((x, m)) 1 ⊕ D κ (x ), • ∀i k + 2, mi = m i-k-2 .
This latter is such that:

G •k+1 g ((x, m)) 1 = G g (G •k g ((x, m)) 1 , mk+1 ) 1 = G g (G •k g ((x, m)) 1 , G •k g ((x, m)) 1 ⊕ D κ (x )) 1 = E κ G •k g ((x, m)) 1 ⊕ G •k g ((x, m)) 1 ⊕ D κ (x ) = E κ G •k g ((x, m)) 1 ⊕ G •k g ((x, m)) 1 ⊕ D κ (x ) = E κ (D κ (x )) =x and G •k+1 g (x, m) 2 = m .
This shows that (x , m ) has been reached starting from B. This fact concludes the proof of the lemma and of the proposition claimed previously.

V. TOPOLOGICAL ENTROPY

Another important tool to measure the chaotic behavior of a dynamical system is the topological entropy, which is defined only for compact topological spaces.

Before studying the entropy of CBC mode of operation, Then, (X , d) is compact must be checked.

A. Compactness study

In this section, (X , d) is a compact topological space will be proven, in order to study its topological entropy later. Firstly, as (X , d) is a metric space, it is separated. It is however possible to give a direct proof of this result: Proposition 6 (X , d) is a separated space.

PROOF Let (x, w) = (x, ŵ) two points of X .

1) If x = x, then the intersection between the two balls B (x, w), 1 2 and B (x, ŵ), 1 2 is empty. 2) Else, it exists k ∈ N such that w k = ŵk , then the balls B (x, w), 10 -(k+1) and B (x, ŵ), 10 -(k+1) can be chosen. Now, the compactness of the metric space (X , d) by using the sequential characterization of compactness will be proven.

Proposition 7 (X , d) is a compact space. PROOF Let X = ((x n , m n )) n∈N be a sequence of X .
There is at least one Boolean vector that appears an infinite number of times in the first components of this sequence, as B N is finite. Let x the lowest of them and I the (infinite) subsequence of X constituted by all the block messages having their first component equal to x.

The first block messages (w n ) 0 of the sequences w n ∈ 0, 2 N -1 N (that are the second components of each couple in the infinite sequence I 0 ) all belong in the finite set 0, 2 N -1 , and so at least one word of this finite set appears an infinite number of times in ((w n ) 0 ) n∈N . Let ω 0 ∈ 0, 2 N -1 be the lowest value occurring an infinite number of times in I, and n 0 the index of its first occurrence, such that x n0 = x, (w n0 ) 0 = ω 0 .

Similarly, the subsequence I 1 of X constituted by the block messages (x n , w n ) such that x n = x and (w n ) 0 = ω 0 is infinite, while all the (w n ) 1 belong in 0, 2 N -1 . So at least one element of 0, 2 N -1 appears an infinite number of times in the second block messages of the second components (w n ) 1 of I 1 . Let ω 1 be the lowest value in 0, 2 N -1 occurring an infinite number of times at this position, and n 1 the index in X of its first occurrence.

A subsequence I 2 = (x n , w n ) of X can be defined again such that ∀n, x n = x, (w n ) 0 = ω 0 , and (w n ) 1 = ω 1 , and a similar argument leads to the definition of ω 2 , the lowest value in 0, 2 N -1 appearing an infinite number of times in the third block messages of the sequences w n ∈ 0, 2 N -1 N of I 3 . This process can be continued infinitely.

Finally the point l = x, (w n k ) k of X is defined; the subsequence (x n k , w n k ) of X converges to l. As for all sequences in X a subsequence that converges in X can be extracted. So, the compactness of X can be concluded.

B. Topological entropy

Let (X, d) be a compact metric space and f : X → X be a continuous map. For each natural number n, a new metric d n is defined on X by

d n (x, y) = max{d(f •i (x), f •i (y)) : 0 ≤ i < n}.
Given any ε > 0 and n 1, two points of X are ε-close with respect to this new metric if their first n iterates are ε-close (according to d).

This metric allows one to distinguish in a neighborhood of an orbit the points that move away from each other during the iteration from the points that travel together. A subset E of X is said to be (n, ε)-separated if each pair of distinct points of E is at least ε apart in the metric d n .

Definition 8 Let H(n, ε) be the maximum cardinality of a (n, ε)-separated set, the topological entropy of the map f is defined by (see e.g., [START_REF] Adler | Topological entropy[END_REF] or [START_REF] Bowen | Entropy for group endomorphisms and homogeneous spaces[END_REF])

h(f ) = lim

VI. CONCLUSION AND FUTURE WORK

In this article, the topological study of the CBC mode of operation has been deepened by evaluating its property of topological mixing. Additionally, a quantitative evaluation has been performed by measuring the level of its topological entropy.

In future work, we will be interesting to show how this characterization of disorder can improve the encryption mode in terms of diffusion and confusion. We will then more specifically focus on the hardware implementation of such modes, and we will check whether adding some properties of disorder will help to face side channel attacks.

  Fig. 1. CBC mode of operation

Theorem 1 n log e n 2 =

 12 Entropy of (X , G g ) is infinite.PROOF Let x, x ∈ B N such that ∃i 0 ∈ 1, N , x i0 = xi0 . Then, ∀w, w ∈ S N , d((x, w); (x, w)) 1 But the cardinal c of S N is infinite, then ∀n ∈ N, c > e n 2. So for all n ∈ N, the maximal number H(n, 1) of (n, 1)-separated points is greater than or equal to e n 2 , and thenh top (G g , 1) = lim 1 n log (H(n, 1)) > lim 1 lim (n) = +∞. But h top (G g , ε) is an increasing function when ε is decreasing, then h top (G g ) = lim ε→0 h top (G g , ε) > h top (G g , 1) = +∞,which concludes the evaluation of the topological entropy of G g .

• The integral part between two points X, Y of the topological space X corresponds to the number of


Abdallah Makhoul

Hassan Harb

Data reduction in sensor networks: performance evaluation in a real environment

Keywords: Wireless sensor networks, data reduction, adaptation and aggregation, Kruskal-Wallis test, telosB mote, real experiments

Data reduction is an effective technique for energy saving in wireless sensor networks. It consists on reducing sensing and transmitting data while conserving a high quality of collected information. In this paper we propose an online data reduction model based on Kruskal-Wallis test that allows sensor nodes to adapt their sensing rates based on the data variance. Then, we propose a local aggregation algorithm to reduce further the data set size before sending to the sink. Experimentation on real telosB sensor network testbed shows the effectiveness of our approach in reducing the size of data transmitted over the network and thus saving energy.

I. INTRODUCTION

W IRELESS sensor networks (WSN) present a low cost solution to enhance our lives. Their main advantages are fast, easy deployment and low maintenance cost [START_REF] Akyildiz | Wireless sensor networks: a survey[END_REF]. Indeed, data reduction is one of the most efficient ways to reduce energy consumption in WSNs. It consists on reducing the amount of data sensed and transmitted to the sink.

In the literature, one can find various data reduction approaches based on in-network processing, data compression or data prediction methods. Adaptive filtering techniques were proposed in [START_REF] Ali | Approximate replication of data using adaptive filters in wireless sensor networks[END_REF], [START_REF] Jain | Adaptive stream resource management using kalman filters[END_REF], [START_REF] Abdolee | Diffusion lms algorithms for sensor networks over non-ideal inter-sensor wireless channels[END_REF], [START_REF] Stern | Processing continuous join queries in sensor networks: A filtering approach[END_REF]. They are based mainly on algorithms like least mean square [START_REF] Jain | Adaptive stream resource management using kalman filters[END_REF], [START_REF] Stern | Processing continuous join queries in sensor networks: A filtering approach[END_REF] and Kalman Filter [START_REF] Ali | Approximate replication of data using adaptive filters in wireless sensor networks[END_REF], [START_REF] Abdolee | Diffusion lms algorithms for sensor networks over non-ideal inter-sensor wireless channels[END_REF]. These works focus on coordinating adaptive filters at the sensor node and prediction techniques at the sink. Each node stops send data when the error between the filter input and the filter output is within a specified threshold. Then, the filters are ready to predict the future data. Other data prediction approaches are also studied to conserve energy in WSN [START_REF] Chu | Approximate data collection in sensor networks using probabilistic models[END_REF], [START_REF] Samarah | Vector-based data prediction model for wireless sensor networks[END_REF], [START_REF] Tulone | Paq: Time series forecasting for approximate query answering in sensor networks[END_REF], [START_REF] Wang | Sensor data modeling and validating for wireless soil sensor network[END_REF], [START_REF] Sinha | Prediction models for energy efficient data aggregation in wireless sensor network[END_REF]. It means to predict future information with the use of various algorithms and prevent transmitting the raw data. Stochastic approaches [START_REF] Chu | Approximate data collection in sensor networks using probabilistic models[END_REF], time-Series forecasting [START_REF] Tulone | Paq: Time series forecasting for approximate query answering in sensor networks[END_REF] and heuristics and algorithms [START_REF] Wang | Sensor data modeling and validating for wireless soil sensor network[END_REF], [START_REF] Sinha | Prediction models for energy efficient data aggregation in wireless sensor network[END_REF] are used. Moreover data compression can be applied in sensor networks. It reduces the size of data transmitted in the network by involving encoding at nodes and decoding at the sink [START_REF] Liang | An efficient and robust data compression algorithm in wireless sensor networks[END_REF], [START_REF] Zimos | Multiterminal source coding with copula regression for wireless sensor networks gathering diverse data[END_REF], [START_REF] He | Compressive data gathering with low-rank constraints for wireless sensor networks[END_REF]. Although these approaches predict sensed values and allow efficient data reduction, however they present several disadvantages. They are almost complex, sometimes they generate communication A. Makhoul is with FEMTO-ST Institute/CNRS, the DISC department, Univ. Bourgogne Franche-Comté, Belfort, France, e-mail: abdallah.makhoul@univ-fcomte.fr.

H. Harb is with the department of Computer Science, American University of Culture and Education (AUCE), Nabatiyeh/Tyre, Lebanon, e-mail: hassanharb@auce.edu.lb. overhead, and the sink may need some transmissions to detect failures.

In this paper, we present two phases data reduction method. It is less complex and suitable for limited resources sensor nodes. The first phase uses Kruskal-Wallis test, instead of Bartlett test used in [START_REF] Makhoul | Residual energy-based adaptive data collection approach for periodic sensor networks[END_REF], and aims to reduce the data acquisition on sensors by adapting their sensing rates to the varying nature of the sensed data. Indeed, Bartlett test has a serious weakness if the data normality assumption is not met. Consequently, in this paper we study the Kruskal-Wallis test which does not impose the data normality assumption. The second phase of our approach consists on aggregating collected data online based on some similarity properties before sending them to the sink. To evaluate our technique, we conducted several experiments on a real environment sensor networks based on telosB nodes. In our experiments, we compared our results to those obtained with Bartlett test [START_REF] Makhoul | Residual energy-based adaptive data collection approach for periodic sensor networks[END_REF] and a data compression method, called S-LEC, proposed recently in [START_REF] Liang | An efficient and robust data compression algorithm in wireless sensor networks[END_REF].

II. DATA ACQUISITION REDUCTION

In this section, we introduce the first phase, i.e. data acquisition phase, which is based on the Kruskal-Wallis test. The idea here is to apply the Kruskal-Wallis test in order to verify if there is high variation in the collected measurements. In the affirmative case the sensing rate must be at its maximum else we adapt the sampling rate according to the variation and to the situation risk as explained in [START_REF] Makhoul | Residual energy-based adaptive data collection approach for periodic sensor networks[END_REF].

Subsequently, the Kruskal-Wallis H test is a rank-based nonparametric test that can be used to determine if there are statistically significant differences between two or more groups of data [START_REF] Mckight | Kruskal-wallis test[END_REF]. We propose to present the Kruskal-Wallis test via the following example.

A. Illustrative Example

Consider the readings taken in three periods p In order to apply the Kruskal-Wallis test, we must first order the readings in all periods by increasing order of their values. The order of each measure represents its rank, denoted r where r ∈ [1, N ]. Then, we assign for each measure its rank as shown in Table I. Subsequently, in case where a measure is repeated several times, called "tied", the mean value of their rank is assigned to the measure. Then, consider the following two hypothesis:

• H 0 : the three probabilities distributions are the same.

• H 1 : the three probabilities distributions are not the same with some desired false-rejection probability (risk α). Then, the Kruskal-Wallis test statistic is given by:

H = 12 N × (N + 1)
r 2 i n i -3 × (N + 1) (1)
where:

• N is the total number of measures in all periods.

• n i is the number of measures during a period.

• r i is the rank of each period.

1 Thus, based on Equation 1, H is calculated as follows:

H = 12 15 × (15 + 1) 43 2 6 + 44 2 5 + 33 2 4 -3 × (15 + 1) = 0.3805
Hence, for α = 0.05 we have H t = 5.991, we have H < H t . Thus, the null hypothesis (H 0 ) is accepted and the sensor sampling rate should be adapted.

Algorithm 1 describes the adaptive sampling rate algorithm at the sensor node based on the Kruskal-Wallis test. For each round, every node decides to increase or decrease its sampling rate according to the difference between its collected measures and the application risk level. First, the node computes the rank for each measure (lines 7-8). Then, it uses the Behavior function [START_REF] Makhoul | Residual energy-based adaptive data collection approach for periodic sensor networks[END_REF] to adapt its sensing rate only if the calculated difference between measures is less than the Kruskal-Wallis threshold (lines 10-14). As explained in [START_REF] Makhoul | Residual energy-based adaptive data collection approach for periodic sensor networks[END_REF], we define the application risk level and we express this level by a quantitative variable R which can take values between 0 and 1 representing the low and the high risk level respectively.

III. SENSORS DATA AGGREGATION A. Definitions and Notations

In periodic applications, each period p is divided into T equal time slots where, at each slot, a sensor S i captures a new reading r ij , then, it forms a vector of readings during the period p as follows: R i = r i1 , r i2 , . . . , r iT .

Mostly, R i may contain redundant (or very similar) readings, especially when the monitored condition varies slowly or when the slots are short. In order to eliminate similar values from the vector R i , we define Similar function as follows:

Definition 1: Similar function. We define the Similar function between two readings as:

Similar(r ij , r i k ) =    1 if |r ij -r i k | = 0, 1 if |r ij -r i k | ≤ δ, 0 otherwise.   
1 the numbers'12' and '3' are fixed by the test and do not depend on the number of periods.

Algorithm 1 Adaptive Sensing Rate Algorithm.

Require: p (1 round = p periods), τ : period size, R: application criticality, α: false-rejection probability. Ensure: S t (instantaneous sampling speed).

1: S t ← τ measures/period 2: while Energy > 0 do 3:

for i = 1 → p do 4:
takes measures at S t speed 5:

end for 6:

for each round do sort measures by increasing order of their values. send (M i ) 18: end while where r ij and r i k ∈ R i and δ is a threshold determined by the application. Furthermore, two readings are considered similar if and only if their Similar function is equal to 1. This means that two measures are considered redundant if they are equals or similar.

Then, we define the weight of a reading as follows: Definition 2: Reading's weight, wgt(r ij ). The weight of a reading r ij is defined as the number of similar readings (according to Similar function) to r ij in the same vector R i .

B. Aggregation Phase Algorithm

Algorithm 2 presents the aggregation process which is running by the sensors themselves at each period. For each new captured reading, S i searches for similarities of the new taken reading. If a similar reading is found, the new one is deleted and the corresponding weight is incremented by 1, else the sensor adds the new reading to the set and initializes its weight to 1. Consequently, S i will possess a set of readings/weights, e.g. R i ={(r i1 , wgt(r i1 )), (r i2 , wgt(r i2 )), . . . , (r i k , wgt(r i k ))} where k ≤ T , which will be sent to the sink.

IV. EXPERIMENTAL RESULTS

To evaluate our proposal, we deployed four Crossbow telosB motes, collecting temperature and humidity measures in our laboratory for about three days (Fig. 1). The proposed approach (Kruskal-Wallis test), the Bartlett test [START_REF] Makhoul | Residual energy-based adaptive data collection approach for periodic sensor networks[END_REF], the S-LEC data compression technique [START_REF] Liang | An efficient and robust data compression algorithm in wireless sensor networks[END_REF] and the naïve approach 2 , e.g. without adapting data collection, are implemented on motes Algorithm 2 Local Aggregation Algorithm.

Require: new reading r i , δ: similarity threshold. Ensure: set of readings with their weights: R i .

1: for each existing reading r j ∈ R i do 2:

if Similar(r i , r j ) = 1 then 3:

wgt(r j ) ← wgt(r j ) + 1 wgt(r i ) ← 1 7:

R i ← R i ∪ {(r i , wgt(r i ))} 8:
end if 9: end for Ids=1, 2, 3 and 4 respectively. Every 30 seconds, each mote collects new measure and sends to a fifth mote, called collector with Id=0, connected to a laptop machine. The collector only relays the data received from the motes to the laptop machine which acts as the sink node. Then, we implemented a Java application on the laptop machine in order to perform daily statistics over the data sent from the motes. In Table II A. Data collected reduction size study Fig. 2 shows a comparison between the number of temperature and humidity measures sent daily by the motes. The obtained results show that Kruskal-Wallis test can reduce up to 19% and 26% of temperature and humidity measures sent with S-LEC technique. Whilst, comparison between Kruskal-Wallis and Bartlett tests shows that for a small value of α, e.g. in days 1 and 2, the mote sends fewer number of measures (temperature and humidity) using Bartlett test, otherwise, data collected are more reduced using Kruskal-Wallis test. This means that Kruskal-Wallis is more flexible than Bartlett to reject the null hypothesis when the false-rejection probability increases.

B. Adapting sensor sensing rate

Fig. 3 shows how the temperature and humidity sensors are able to adapt their sampling rate using the Kruskal-Wallis test. The obtained results show that the sampling rate of each sensor is dynamically adapted after each round. This confirms the decrease in the volume of the collected data compared to the naïve approach, e.g. S t = 50 measures/period. We can also notice that the sampling rate of the temperature sensor has been more adapted compared to the humidity sensor.

Fig. 5 shows an illustrative example for which data are collected by the temperature sensor using Kruskal-Wallis test. Compared to naïve collection, the results show that our adapting algorithm allows each sensor to collect most important measures in the round after adapting its sampling rate.

C. Energy consumption study

In nesC programming, there is no model provided by tinyOS in order to measure the energy consumed in telosB [16]. In our experiments, we used the same radio model used in [START_REF] Madden | Intel berkeley research lab[END_REF] as the most used model to evaluate the energy consumption in WSNs. Since the energy consumption is highly related to the amount of data collected/sent (see Fig. 2), the obtained results show that: a) Kruskal-Wallis test conserves the energy of the mote up to 44% compared to S-LEC. b) for α equals to 0.01 and 0.025, the Bartlett test consumes about 12% less energy compared to the Kruskal-Wallis test in the first two days. c) For α = 0.05 in the third day, Kruskal-Wallis test reduces the energy consumption in the mote about 6% compared to Bartlett test.

D. Data accuracy/integrity

Fig. 6 shows the percentage of data loss using the three approaches. The percentage of data loss is calculated by searching the measures collected in the naïve approach (mote Id=1) and do not collected (nor similar values) by the other motes. The obtained results are dependent on the number of data sent to the sink (see Fig. 2); thus, more data are sent less measures are lost. Indeed, we observe that the percentage of data loss using kruskal-Wallis test does not accessed 3.1% whereas it arrives to 6.4% and 4.3% using S-LEC and Bartlett test. Therefore, we consider that Kruskal-Wallis test is more efficient in terms of saving data integrity compared to other approaches.

E. Further Discussion

In this section, we make more discussion about the three compared approaches then, we give some directions as to which approach is more suitable for a desired application. Regarding the energy preserving viewpoint, Kruskal-Wallis and Bartlett tests significantly reduce the energy consumption in motes, compared to S-LEC approach. Furthermore, comparison both tests shows that Bartlett conserves more energy when α is small (e.g. ≤ 0.025) whereas, the energy is more preserved using Kruskal-Wallis for a higher value of α (e.g. > 0.025). Therefore, for applications where we need to conserve the energy of the network, the decision makers should be used the suitable test depending on the chosen value of α. Regarding the data accuracy viewpoint, Kruskal-Wallis test can save the integrity of the collected data more than those in other approaches in almost cases. Hence, it is suitable to use the Kruskal-Wallis test for the application that does not permit flexibility regarding data accuracy. V. CONCLUSIONS In this paper, we proposed energy-efficient data reduction and aggregation techniques dedicated to wireless sensor networks. They were targeted to minimize the amount of data retrieved/communicated by the network without loss in fidelity. We studied a data collection model based on a Kruskal-Wallis test that allows each sensor node to adapt its sensing rate to the changing of the monitored condition. Furthermore, we added a second level of local data aggregation in order to reduce data transmitted in the network and save energy. We showed via real experimentation and telosB testbed that our approach can be effectively used to increase the sensor network lifetime, while still keeping the quality of the collected data high.

As a future work, we seek to extend our adaptive sensing technique in order to take into account the correlation between neighboring sensor nodes.

Fig. 1 .Fig. 2 .Fig. 3 .Fig. 4 .Fig. 5 .Fig. 6 .

 123456 Fig. 1. Sample of TelosB nodes used in the experimentation.

  1 , p 2 and p 3 as shown in the following table.

		p1	r1	p2	r2	p3	r3	Total
		8.2	1	10.2	7	13.5	12
		10.3	8	9.1	4→ 4	8.4	2
		9.1	3→ 4	13.9	14	9.6	6
		12.6	10	14.5	15	13.8	13
		11.4	9	9.1	5→ 4
		13.2	11
	ni	6		5		4		N = 15
	T otal		43		44		33

  S t ← BV (H, H t , R, τ ) (BV behavior function).

	8:	compute the rank of each measure
	9:	find H t
	10:	if H ≤ H t then
	11:
	12:	else
	13:	S t ← τ measures/period
	14:	end if
	15:	end for
	16:	M

i ← local Aggregation(M i , δ) // M i is the set of measures collected at the current period 17:

  , we show the metrics used in our experiments.

	Parameter	Description	Value
	round	round size	2 periods
	τ	period size	50 measures
	δ	similarity threshold	0.05 for temperature
			0.01 for humidity
	R	application criticality	0.6
	α	false-rejection probability	0.01 in day 1
			0.025 in day 2
			0.05 in day 3

TABLE II .

 II EXPERIMENTAL ENVIRONMENT.

naïve approach has been implemented in order to calculate the data loss measures in the two compared tests.


Stéphane Chrétien

Sébastien Darses
email: sebastien.darses@univ-amu.fr

Christophe Guyeux
email: christophe.guyeux@univ-fcomte.fr

Paul Clarkson
email: paul.clarkson@npl.co.uk

Stephane Chretien
email: stephane.chretien@npl.co.uk

On the pinning controllability of complex networks using perturbation theory of extreme singular values. Application to synchronisation in power grids

Keywords: 1991 Mathematics Subject Classification. Primary: 93Bxx, Secondary: 15A18 Pinned control, eigenvalue perturbation, sum of rank one matrices, singular value perturbation, control theory

, where it was reformulated as a global asymptotic stability problem. The goal of this short note is to refine the analysis proposed in [16] using recent results in singular value perturbation theory.

Introduction

In recent years, extensive efforts have been devoted to the control of complex dynamical networks [13,[START_REF] Slotine | Applied nonlinear control[END_REF]. Real networks for which such problems occur are of paramount importance in the natural sciences and engineering [START_REF] Porfiri | Criteria for global pinning-controllability of complex networks[END_REF]. In particular, the examples of contemporaneous beats of the heart cells [25] or the synchronous behaviors of the cells of the suprachiasmatic nucleus in the brain [17] have been considered in the literature. Social networks are also fascinating examples where the formation of mass opinions and the emergence of collective behaviors are frequently observed. One major issue is that real world networks usually consist of a very large number of nodes and links which makes it impossible to apply control actions to all nodes.

Pinning control is a new way to address such control problems by placing local feedback injections on a small fraction of the nodes. Controllability of such systems has been addressed in [START_REF] Porfiri | Criteria for global pinning-controllability of complex networks[END_REF], where it was reformulated as a global asymptotic stability problem. The goal of this short note is to refine the analysis proposed in [START_REF] Porfiri | Criteria for global pinning-controllability of complex networks[END_REF] using recent results in singular value perturbation theory.

1.1. Notations. The Kronecker symbol is denoted by δ i,j , i.e., δ i,j = 1 if i = j and is equal to zero otherwise. We denote by x 2 the euclidian norm of a vector x and by A the associated operator norm (spectral norm) of a matrix A. Its transposition is denoted by A t .

For any symmetric matrix B ∈ R d×d we will denote its eigenvalues by

λ 1 (B) ≥ • • • ≥ λ d (B).
The largest eigenvalue will sometimes also be denoted by λ max (B) and the smallest by λ min (B). The smallest nonzero eigenvalue of a positive semi-definite matrix B will be denoted by λ min>0 (B). Finally, deg i is the degree of a vertex i in a given graph, that is, its number of edges.

1.2. The model. One considers a set of N n-dimensional oscillators governed by a system of nonlinear differential equations. Moreover, we assume that each oscillator is coupled with a restricted set of other oscillators. This coupling relationship can be efficiently described using a graph where the vertices are indexed by the oscillators and there is an edge between two oscillators if they are coupled. The overall dynamical system is given by the following set of differential equations

x i (t) = f (x i (t)) -σB N j=1 l ij x j (t) + u i (t), t ≥ t 0 , (1.1) i = 1, . . . , N , where x i (t) ∈ R n is the state of the i th oscillator, σ > 0, B ∈ R n×n , f : R → R describes the dynamics of each oscillator, L = (l ij ) i,j=1,...,N
is the graph Laplacian of the underlying graph, and u i (t), i = 1, . . . , N are the controls. For the system to be well defined, we have to specify some initial conditions x i (t 0 ) = x i0 for i = 1, . . . , N .

1.3. The control problem. Assume that we have a reference trajectory s(t), t ≥ t 0 satisfying the differential equation

s (t) = f (s(t)).
Our goal is to control the system using a limited number of nodes. The selected nodes are called the "pinned nodes". For this purpose, we use a linear feedback law of the form

u i (t) = p i Ke i (t),
where e i (t) = s(t) -x i (t), K is a feedback gain matrix, and where

p i = 1 if node i is pinned 0 otherwise.
Let P denote the diagonal matrix with diagonal (p 1 , . . . , p N ).

1.4. Controllability. In [START_REF] Porfiri | Criteria for global pinning-controllability of complex networks[END_REF], the authors propose a definition for (global pinning-) controllability (based on Lyapunov stability criteria):

Definition 1.1. We say that the system (1.1) is controllable if the error dynamical system e := (e i (t)) 1≤i≤N is Lyapunov stable around the origin, i.e., there exists a positive definite function V such that d dt V (e(t)) < 0 when e(0) = 0.

The following result, [START_REF] Porfiri | Criteria for global pinning-controllability of complex networks[END_REF]Corollary 5], provides a sufficient condition for a system to be controllable: [START_REF] Porfiri | Criteria for global pinning-controllability of complex networks[END_REF]). Assume that f is such that there exists a bounded matrix F ξ, ξ , whose coefficients depend on ξ and ξ, which satisfies

Proposition 1.2 ([
F ξ, ξ ξ -ξ = f (ξ) -f ( ξ), ξ, ξ ∈ R n . (1.2)
Let Q ∈ R n×n be a positive definite matrix such that

QK + K t Q t = κ QB + B t Q t QB + B t Q t 0 and 1 2 λ N (σL + κP ) λ n QB + B t Q t > sup ξ, ξ F ξ, ξ Q . (1.3)
Then the system is controllable.

Many systems of interest satisfy the constraint specified in Proposition 1.2, see [START_REF] Jiang | A simple global synchronization criterion for coupled chaotic systems[END_REF]. This proposition is very useful for node selection via the matrix P . Indeed, assume that Q is selected, then one may try to maximise λ N (σL + κP ) as a function of P , under the constraint that no more than r nodes can be pinned. This is a combinatorial problem that can be relaxed using semi-definite programming or various heuristics [START_REF] Ghosh | Growing well connected graphs[END_REF].

1.5. Goal of the paper. Our goal in the present note is to propose an easy controllability condition refining [16, Corollary 7], based on the algebraic connectivity of the graph, the number of pinned nodes, the coupling strength and the feedback gain. Our approach is based on perturbation theory of the extreme singular values of a matrix after appending a column. The basic results of this theory are given in the appendix.

Main result

2.1.

A simple criterion for controllability. Our main result is the following theorem.

Theorem 2.1. Let Q ∈ R n×n be a positive definite symmetric matrix that satisfies

QK + K t Q t = κ QB + B t Q t QB + B t Q t 0,
and assume that

σλ min>0 (L) λ min QB + B t Q t > sup ξ, ξ F ξ, ξ Q . (2.4) If κ satisfies κ ≥ r i=1 deg i σλ min>0 (L) - 2 F ξ, ξ Q λ min (QB+B t Q t ) + σλ min>0 (L), (2.5)
then the system is controllable.

Proof. We follow the same steps as for the proof of Corollary 7 in [START_REF] Porfiri | Criteria for global pinning-controllability of complex networks[END_REF]. We assume without loss of generality that the first r nodes are the pinned ones. We may write P as

P = r i=1 e i e t i ,
where e i is the i th member of the canonical basis of R N , i.e., e i (j) = δ i,j . We will try to compare λ N (σL + κP ) with λ N (σL) and use Proposition 1.2 to obtain a sufficient condition for controllability based on L, i.e. the topology of the network. For this purpose, let us recall that L can be written as

L = I • I t ,
where I is the incidence matrix of any directed graph obtained from the system's graph by assigning an arbitrary sign to the edges [START_REF] Brouwer | Spectra of graphs[END_REF]. Of course L will not depend on the chosen assignment. Using this factorization of L, we obtain that σL + κ Moreover, λ min>0 (σL + κP ) can be expressed easily as the smallest nonzero eigenvalue of the r th term of a sequence of matrices with shape (A. [START_REF] Porfiri | Criteria for global pinning-controllability of complex networks[END_REF]) for which we can use Theorem A.2 iteratively. Indeed, we have

λ min>0 (σL + κe 1 ) = λ min>0 √ κ e 1 , √ σI t √ κ e 1 , √ σI .
Let us denote by x the vector √ κ e 1 and by X the matrix [ √ σI]. Then, we have that

√ κ e 1 , √ σI t √ κ e 1 , √ σI = x t x x t X X t x X t X .
Therefore, Theorem A.2 gives

λ min>0 σL + κe 1 e t 1 ≥ σλ min>0 (L) - deg 1 (κ -σλ min>0 (L))
, .

where
Since λ min>0 (σL + κe 1 e t 1 ) ≤ λ min>0 (σL), we thus obtain

λ min>0 σL + κe 1 e t 1 + κe 2 e t 2 ≥ λ min>0 (σL + κe 1 e t 1 ) - deg 2 (κ -σλ min>0 (L))
.

We can repeat the same argument r times and obtain

λ min>0 (σL + κP ) ≥ σλ min>0 (L) - r i=1 deg i κ -σλ min>0 (L) . (2.6)
Finally, by Proposition 1.2, we know that the following constraint is sufficient for preserving controllability

λ min>0 σL + κ r i=1 e i e t i ≥ 2 F ξ, ξ Q λ min (QB + B t Q t ) . (2.7)
By (2.6), it is sufficient to guarantee the controllability of our system to impose

σλ min>0 (L) - r i=1 deg i κ -σλ min>0 (L) ≥ 2 F ξ, ξ Q λ min (QB + B t Q t ) .
Then, combining (2.7) with (2.4) implies that

κ ≥ r i=1 deg i σλ min>0 (L) - 2 F ξ, ξ Q λ min (QB+B t Q t ) + σλ min>0 (L)
is a sufficient condition for controllability.

2.2.

Consequence of the formula. What the formula (2.5) tells us is that the nodes of smallest degree should be selected before the others. The importance of considering the smallest degrees first comes as a natural consequence of our formula in a very explicit fashion, and the proof of the formula is slightly more transparent than the argument provided, e.g., in [START_REF] Porfiri | Criteria for global pinning-controllability of complex networks[END_REF].

2.3.

Comparison with previous works. Closely related statements in the context of numerical and theoretical optimization of pinning strategies have already been established in the literature (cf. [START_REF] Porfiri | Experiments on node-to-node pinning control of Chua's circuits[END_REF] and [START_REF] Sorrentino | Effects of the network structural properties on its controllability[END_REF]), while in [START_REF] Delellis | Pinning control of complex networks via edge snapping[END_REF] a more recent work on edge snapping provides interesting insight into optimal evolution of pinning gains. For instance, the idea of using the node degrees for choosing the pinned nodes of the network already appeared in [START_REF] Porfiri | Criteria for global pinning-controllability of complex networks[END_REF] but not as explicitly as in our work, as it is a consequence of a subsequent qualitative discussion. Let us also mention that the importance of considering the degree sequence of the graph also appears in the more recent works [START_REF] Porfiri | Node-to-node pinning control of complex networks[END_REF] but only via the maximum and minimum degree and [START_REF] Sorrentino | Effects of the network structural properties on its controllability[END_REF], where the analysis is done via computational experiments.

Additionally, we used a perturbation theorem that is more precise than Weyl's perturbation bound, and therefore, we end up with tighter inequalities. Finally, κ can be explicitly found with our approach, while it is implicit in [START_REF] Porfiri | Criteria for global pinning-controllability of complex networks[END_REF].

Application to synchronisation in power grids

Our methods can be applied to power network control and synchronization as described in [START_REF] Dörfler | Synchronization in complex networks of phase oscillators: A survey[END_REF][START_REF] Nishikawa | Comparative analysis of existing models for power-grid synchronization[END_REF][START_REF] Motter | Spontaneous synchrony in power-grid networks[END_REF]. We consider an AC power network modelled by a graph G = (V, E), n = |V| where the vertices represent the nodes and the edges represent the transmission lines. The network is also defined by a symmetric admittance matrix Y .

The circuit equations have a solution V i at each node i = 1, . . . , N . The set V of nodes is divided into three different subsets. The set V 1 contains the index set of the load buses, V 2 is the index set of the synchronous generators, and V 3 is the index set of the grid connected direct current (DC) power sources. For lossless networks, the active power flow from node j to node j is given by a jj sin(θ j -θ j ) where θ j is the argument of V j and

a jj = |V j | |V j | I(Y jj ). (3.8)
The system is governed by the following equations. For load buses, we have

D j dθ j dt + P (1) j = - N j =1
a jj sin(θ j -θ j ), j ∈ V 1 , (3.9) where P

(1) j > 0 is the active power drawn at node j (1 in the exponent is a notation to recall that we consider the set V 1 ) and D j > 0 is the damping factor. In the case of synchronous generators, we have

M j d 2 θ j dt 2 + D j dθ j dt = P (m) j - N j =1
a jj sin(θ j -θ j ), j ∈ V 2 , (3.10) where P (m) j > 0 is the mechanical power output (m in exponent is for mechanical ) and M j > 0 is the inertia factor. Finally, for DC power sources V 3 , we have:

D j dθ j dt = P (d) j - N j =1 a jj sin(θ j -θ j ), j ∈ V 3 , (3.11)
where we assumed that each DC source is connected to the AC grid via and DC/AC inverter and that each inverter j ∈ V 3 has a nominal value P (d) j > 0 and is equipped with a droop controller with droop slope 1/D j , following the model of [START_REF] Simpson-Porco | Synchronization and power sharing for droop-controlled inverters in islanded microgrids[END_REF].

After linearisation, we obtain the system

dx j dt = - 0 P (1) j D j + n j =1, j =j 1 D j 0 0 0 a jj x j -   n j =1, j =j 1 D j 0 0 0 a jj   x j , j ∈ V 1 , (3.12
)

dx j dt =   0 P (m) 2 M 2   N j =1, j =j 0 0 a 2,j M 2 0 x j -   N j =1 0 0 a j,j M 2 0 +   0 D 2 M 2 0 0   +   x j , j ∈ V 2 , (3.13)
and

dx j dt = - 0 P (1) j D j + n j =1, j =j 1 D j 0 0 0 a jj x j -   n j =1, j =j 1 D j 0 0 0 a jj   x j , j ∈ V 3 . (3.14)
Denote by P the vector whose j th component is P

j if j ∈ V 1 , P (m) j if j ∈ V 2 , (1)
otherwise. The results of the previous section readily apply to the system described by (3.12), (3.13) and (3.14) if we duplicate the components of P in order to take into account that the state at each bus is two-dimensional.

Figures 1,2, and 3 show the evolution of the phase at each node for the system governed by (3.12), (3.13), and (3.14). The way we obtained them is described in details hereafter. As stated previously, we implemented our method in Python. In particular, we used the NetworkX [START_REF] Aric | Exploring network structure, dynamics, and function using NetworkX[END_REF] library for graph modelling, we used a SciPy [START_REF] Jones | SciPy: Open Source Scientific Tools for Python[END_REF] library module for simulating the power grid's dynamics, and Numpy [START_REF] Stéfan Van Der Walt | The NumPy Array: A Structure for Efficient Numerical Computation[END_REF] to compute with matrices and generate random objects. We generated the network based on 3 main parameters:

• the total number n of nodes (or buses, n = 50 in the simulations shown in the figures), • a threshold T under which the coupling strength is ignored (set to 0 in the figures), and • the number of nodes (ordered by increasing degree) that are pinned (10 here). The vectors P , D, M and K which appear in formulas (3.9), (3.10), and (3.11) are chosen at random in [0,1]. The only components of K considered in the simulations are the ones associated by the pinned nodes. The coupling strength matrix a is random with i.i.d. components in [0,1]. This latter is used as connectivity matrix, except that edges with weight lower than T are removed. The n nodes are then sorted in increasing degree order, and divided into three sets, namely the V 1 , V 2 , and V 3 . The ODE's are then defined and solved numerically using the "integrate" SciPy module. The Figure 1. Evolution of the phase at each bus as a function of time for a random network integrator chosen for these experiments was vode [START_REF] Brown | VODE: A Variable C Coefficient ODE Solver[END_REF], which is a real-valued variable-coefficient ODE solver with fixed-leading-coefficient implementation, with a backward differentiation formulas method. The iteration method of the ODE solver's correction step is "chord iteration" with an internally generated full Jacobian. The maximum order used by the integrator was set to 5, while maximum number of (internally defined) steps allowed during one call to the solver was equal to 500,000. Finally, evolution of the phase at each bus has been plotted using matplotlib [START_REF] Hunter | Matplotlib: A 2D graphics environment[END_REF] library on data produced by the integrated ode object. As depicted in the 3 figures provided below, the method successfully controls the random power grid [START_REF] Abaid | Consensus over numerosity-constrained random networks[END_REF] in this randomly drawn example. We observed that the procedure successfully solved the control problem in most random instances.

Conclusion

In this article, an easy controllability condition refining [16, Corollary 7], based on the algebraic connectivity of the graph, the number of pinned nodes, the coupling strength, and the feedback gain, has been proposed. The result is based on perturbation theory of the extreme singular values of a matrix after appending a column. The method has finally been applied to power network control and synchronization, showing that the proposed results can readily be applied to such systems.

The application considered in the paper, which was an important motivation for our results, still raises important remaining questions. The most difficult is the one of how to handle nonlinear settings of the type arising in power grids. Future work will be devoted to this problem with many expected potential application to other types of networks as well. In particular, the following perturbation bounds of Weyl and Mathias hold:

Corollary A.

  r , . . . , √ κ e 1 , √ σI √ κ e r , . . . , √ κ e 1 , √ σI t .

Figure 2 .Figure 3 . 3 .a 2 η r + η 2 r + 4 a 2 ,

 233222 Figure 2. A second example of evolution of the phase at each bus as a function of time for a random network, with same setup

3 .λ 2 2

 32 r+1 (A) ≥ min(c, λ r ) -a 2 (A.19) λ r+1 (A) ≥ min(c, λ r ) -a |c -λ r | . (A.20)

  deg 1 is the degree of node number 1. Let us now consider λ min>0 (σL + κ e 1 + δ 2 e 2 ). We have that λ

	Let us denote by x the vector	√ κ e 2 and by X the matrix [ √	κ e 1 ,	√	σI].
	Then, we have that
	√	κ e 2 ,	√	κ e 1 ,	√	σI	t √	κ e 2 ,	√	κ e 1 ,	√	σI =	x t x x t X X t x X t X
	and using Theorem A.2 again, we obtain
	λ min>0 σL + κe 1 e t 1 + κe 2 e t 2	≥ λ min>0 (σL + κe 1 e t 1 ) -	deg 2 (κ -λ min>0 (σL + κe 1 e t 1 ))
									min>0	√	κ e 2 ,	√	κ e 1 ,	√	σI	t √ κ e 2 ,	√	κ e 1 ,	√	σI .

min>0 (σL + κ e 1 + δ 2 e 2 ) = λ

Acknowledgement -We thank Maurizio Porfiri for pointing out a missing assumption in Theorem 2.1. The authors also thank the anonymous referees for their constructive and insightful suggestions.

the singular values of the matrix (x, X) obtained by appending the column vector x to the matrix X:

(A1) Consider the matrix

On one hand, one may study in (A1) the eigenvalues of the (n+1)×(n+1) hermitian matrix A, i.e., the matrix X t X augmented with an arrow matrix.

On the other hand, one will deal in (A2) with the eigenvalues of the d × d hermitian matrix A, which may be seen as a rank-one perturbation of XX t . The matrices A and A have the same non-zeros eigenvalues, and in particular λ max (A) = λ max ( A). Moreover, the singular values of the matrix (x, X) are the square-root of the eigenvalues of the matrix A.

Equivalently, the problem of a rank-one perturbation can be rephrased as the one of controlling the perturbation of the singular values of a matrix after appending a column. A.2. A theorem of Li and Li. In this paper, we use a slightly more general framework than (A1), that is the case of a matrix A = c a t a M , (A. [START_REF] Porfiri | Criteria for global pinning-controllability of complex networks[END_REF] where a ∈ R d , c ∈ R and M ∈ R d×d is a symmetric matrix.

The following theorem provides sharp upper bounds for λ max (A), and lower bounds on λ min (A), depending on various information on the submatrix M of A. As discussed above, this problem has close relationships with our problem of appending a column to a given rectangular matrix, because λ 1 ( A) = λ 1 (A).

Theorem A.1 (Li-Li's inequality and a lower bound). Let d be a positive integer and let M ∈ C d×d be an Hermitian matrix, whose eigenvalues are

Let A be given by (A.16). Therefore:

≤ λ 1 (A) -max(c, λ 1 ) ≤ 2 a 2


Suyao Kong
email: suyao.kong@utbm.fr

Mickaël Hilairet
email: mickael.hilairet@univ-fcomte.fr

Robin Roche Advanced

R Roche

Advanced passivity-based control for a

Keywords: Energy Management, Fuel Cell, IDA-PBC, Passivity Based Control, Supercapacitors

come

INTRODUCTİON

Interconnection and Damping Assignment-Passivity Based Control (IDA-PBC) is applied in a fuel cell/super-capacitor hybrid power system. This nonlinear controller uses matrices that are related to the interconnection between the subsystems and the damping of the system to find a desired command that achieves the stability of the whole closed-loop system. M. Hilairet et al. [START_REF] Hilairet | A passivity-based controller for coordination of converters in a fuel cell system[END_REF] applied IDA-PBC for a regular two-converter parallel system with a fuel cell (FC) and supercapacitors (SCs). In their research, some terms of the general control were set to zero in order to obtain regular control strategies.

The novelty of this work is to exploit all the terms of the general non-linear control in order to integrate some components constraints directly in the control, while preserving the stability of the whole closed-loop system.

ADVANCED POWER MANAGEMENT STRATEGY

A two-converter parallel structure is used for the fuel cell/super-capacitors system shown in Fig. 1.

Fig. 1: System model

The main idea of the energy management strategy is to use the FC to supply the load and use the SCs during transients in order to stabilize the DC bus, with a proof of the stability of the whole closed-loop system. Based on the previous strategy described in [START_REF] Hilairet | A passivity-based controller for coordination of converters in a fuel cell system[END_REF], a saturation of the SCs current is added to consider the state-of-charge of SCs by integrating a variable 𝑟 2 . The control law is as follows:

𝑖 𝑠𝑐 * * = 𝑖 𝑠𝑐 * + 𝑟 2 𝐶 𝑠𝑐 2 𝑣 ̃𝑠𝑐 (1)
𝑖 𝑓𝑐 * * = 𝑖 𝑓𝑐 * -𝑟 2 𝑣 𝑠𝑐 𝑣 𝑓𝑐 𝐶 𝑠𝑐 2 𝑣 ̃𝑠𝑐 (2)
where 𝑖 𝑠𝑐 * , 𝑖 𝑓𝑐 * are control inputs obtained in [START_REF] Hilairet | A passivity-based controller for coordination of converters in a fuel cell system[END_REF], and 𝑣 ̃𝑠𝑐 is the error between the voltage of super-capacitors and its equilibrium point.

SİMULATİON RESULTS

Fig. 2 shows a comparison between the original control (in blue) proposed in [START_REF] Hilairet | A passivity-based controller for coordination of converters in a fuel cell system[END_REF] and the new one (in red). Here, the voltage of SCs is limited between 20.5 V and 21.5 V for the purpose of the test. It follows that the SCs current goes to 0 when the voltage exceeds the limitation around 20 and 30s.

In order to counteract the injection of current in the DC bus, the FC current is naturally increased, and consequently preserves the stability of the whole closed-loop system.

Fig. 2 :

 2 Fig. 2: Simulation results