Isabelle Oliveira
email: ioliveira@isc.cnrs.fr

Sabine Ploux
email: sploux@isc.cnrs.fr

Keywords: métaphore, traitement automatique des langues, corpus

Vers une méthode de détection et de traitement automatique de la métaphore

Résumé

Nous présentons un modèle pour le repérage automatique de la métaphore et pour l'aide à son analyse. Les résultats actuels de ce modèle sont très prometteurs. La présentation du modèle est suivie d'une analyse comparative de la métaphore en langue de spécialité (ici le domaine du droit) et en littérature. Nous projetons de confronter ces deux approches afin à la fois de valider systématiquement le système d'analyse automatique et de donner une base objective à l'étude linguistique de la métaphore.

Introduction

La métaphore terminologique, processus linguistique fondé sur le rapprochement mental qu'opère l'utilisateur d'une langue entre deux objets ou phénomènes par association d'idées, est une arme précieuse dans l'arsenal linguistique du juriste, mais aussi du littéraire. Cependant, en langue de spécialité on cherche à se démarquer de la métaphore qui renvoie souvent à la polysémie et à l'ambiguïté. En effet, cette figure de style consistant en un transport de sens risque à tout moment de rendre inopérant l'effort d'univocité et de rigueur, mais la métaphore n'est pourtant pas absente des langues de spécialité. Il faut alors se demander si cette présence est accidentelle, contingente, ou si elle est incontournable et signifiante. Afin d'étudier ces questions, nous avons cherché dans un premier temps à nous donner les moyens de procéder à une analyse systématique effectuée sur des corpus de textes. Dans la première partie de l'article nous présenterons, l'ensemble des corpus constitués et les méthodes de repérage et d'analyse des métaphores qui en sont extraites. Dans la seconde partie, nous présentons une analyse linguistique comparative de la métaphore littéraire et de la métaphore dans une langue de spécialité : le droit. Notre travail est prospectif et nous projetons dans l'avenir de rapprocher et de faire interagir ces deux facettes de l'analyse afin à la fois d'affiner la démarche automatique et de nourrir l'étude linguistique de résultats obtenus non plus « à la main » mais par un traitement systématique et automatique.

2 Une approche fondée sur l'analyse de corpus et le traitement automatique des langues

Description des corpus

Les corpus que nous avons compilés ont été regroupés suivant leur genre : littéraire, journalistique et de droit. Du point de vue de l'étude de la métaphore, ces trois genres présentent des caractéristiques différentes qui, comme nous le montrerons dans la seconde partie de l'article, justifient leur choix. Chaque corpus contient des textes en français et des textes en portugais. Nous nous sommes efforcées, quand cela était possible, de constituer des corpus parallèles et sinon comparables. Ce choix permettra par la suite une étude translingue du phénomène métaphorique. Corpus journalistique Le corpus journalistique contient 5 années du journal Le Monde (1997Monde ( -2001)). Nous procédons actuellement à la constitution d'un corpus comparable à partir de quotidiens portugais.

Corpus littéraire Français-Portugais

Corpus de droit

Le corpus parallèle de droit contient les traités communautaires (Constitution européenne, traités d'Amsterdam, de Lisbonne, de Maastricht, de Madrid, de Rome et le traité constituant la Communauté européenne de l'énergie atomique (CEEA)) dans les deux langues française et portugaise.

Vers une détection automatique de la métaphore

Le repérage automatique ou semi-automatique des métaphores constitue sans doute à la fois un défi majeur pour le traitement automatique des langues et une perspective utile pour le travail linguistique. Dans ce projet, nous avons cherché un critère utile à ce repérage. Ce critère est fondé sur la définition même de la métaphore à savoir l'association entre des domaines source et cible hétérogènes et repose sur une mesure de la distance sémantique entre les domaines en interaction. Ainsi, dans la métaphore juridique fruit civil, les domaines hétérogènes de la botanique et du droit sont liés et associés pour produire le sens voulu. A contrario, dans l'expression fruits et légumes, les deux termes (fruit et légume) ont des contextes d'emploi très souvent similaires. Dans ce travail, nous proposons d'utiliser le modèle ACOM (Automatic Contexonym Organizing Model) [START_REF] Ji | Lexical knowledge representation with contexonyms[END_REF] des Atlas sémantiques (AS) pour calculer une distance entre les contextes d'emploi les plus génériques des termes d'une expression (métaphorique ou non) et nous utilisons le résultat de ce calcul comme critère de détection d'une métaphore.

Le modèle ACOM, dont on peut consulter les résultats à l'adresse http://dico.isc.cnrs.fr, fournit pour chaque mot une liste de mots associés qui représentent et caractérisent son emploi. Il est fondé sur une analyse systématique de corpus et a été testé pour la valeur cognitive de ses résultats. En effet, les mots proposés par ACOM en réponse à une requête sont souvent identiques aux mots donnés par des sujets humains lors d'une tâche d'association [START_REF] Ji | Testing the Cognitive Relevance of a Geometric Model on a Word-Association Task: A Comparison of Humans, ACOM, and LSA[END_REF]. Ainsi les premiers mots associés au mot poire et donnés par le modèle sont abricot, chou-fleur, couper, morceaux, doigts, fromage, fruits, pêche... On pourra se reporter à l'article mentionné ci-dessus pour plus de précisions concernant le modèle et sa mise en oeuvre. Dans la présente étude, nous l'avons utilisé de la façon suivante. Nous avons voulu mesurer le recouvrement thématique entre les termes d'une expression métaphorique par un calcul du recouvrement entre les mots associés à chacun d'eux. L'hypothèse sous-jacente consiste à prédire que dans une métaphore le recouvrement sera réduit car les domaines source et cible sont par définition différents alors que dans une expression non métaphorique ce recouvrement devrait être plus important car les mots associés sont au moins partiellement partagés.

Description de l'expérience réalisée sur corpus

Afin de tester cette hypothèse nous avons constitué des bases de syntagmes nominaux et verbaux extraits du corpus littéraire français et du corpus journalistique français pour en extraire des expressions porteuses de métaphores et des ensembles contrôle d'expressions non métaphoriques.

Constitution des bases de syntagmes à partir des corpus La méthode de constitution de ces bases est décrite dans [START_REF] Ben | Vers un modèle pour la représentation de la sémantique des expressions dans le paradigme géométrique des Atlas sémantiques[END_REF]

Vers une aide à l'analyse automatique de la métaphore

Outre le repérage, nous cherchons une méthode automatique pour l'aide à la détermination des éléments qui participent au processus métaphorique. Un premier exemple de cette démarche a été présenté en 2008 (Ploux et [START_REF] Oliveira | Pour un traitement automatique de la métaphore[END_REF]. En effet, dans l'analyse d'une métaphore, la mise en évidence des traits sémantiques opérants est un travail d'experts dont le résultat varie suivant les auteurs. Il est alors intéressant de chercher une méthode objective qui ne repose pas uniquement sur des interprétations, interprétations qui ne sont ni toujours congruentes ni toujours exhaustives. Notre ambition est alors de prendre support sur une analyse systématique de corpus pour l'élaboration d'un modèle capable de mettre en évidence automatiquement des traits sémantiques à l'oeuvre dans une métaphore. La proposition encore embryonnaire que nous faisons ici repose comme précédemment sur le modèle géométrique ACOM. Ji et ses collègues [START_REF] Ji | Lexical knowledge representation with contexonyms[END_REF] ont montré que ce modèle permet, outre l'utilisation qui en a été décrite au paragraphe précédent, de calculer, pour une liste de mots donnée, des termes qui représentent les contextes d'emploi qui lient les mots de la liste. À titre d'exemples, nous donnons ici pour quelques métaphores (1) les premiers termes associés à chacun des termes de l'expression métaphorique pris isolément ainsi que (2) les premiers termes associés à l'ensemble de la métaphore. On voit sur ces exemples se dégager quelques traits sémantiques qui participent à la cohésion métaphorique. La métaphore est aussi le lien d'une thématique de transfert qui dynamise les représentations construites par les textes. Déplacements corporels, métamorphoses, animation des inanimés de façon pittoresque, sont autant de thèmes récurrents qui constituent l'univers littéraire. La métaphore en littérature lie le sensible et l'intelligible. En voyant ou en lisant, on comprend. Elle n'est rien d'autre que ce miroir qui nous anime par le mouvement ayant pour caractéristique principale d'être créative. Il y aurait ainsi une dynamique propre à la métaphore littéraire qui inviterait à poursuivre l'analogie ou à la critiquer, c'est-à-dire à penser puisque nous voyons ici en quel sens l'on peut dire que pensée et langage ne sont pas deux temps séparés mais les deux faces d'un même processus de recherche.

 Une foule tiède

Spécificité de la langue juridique

La langue juridique s'associe plus à la logique, qui vise le repérage d'entités absolument stables, non affectées par l'incessante ambigüité, ambivalence, qui apparaissent ou disparaissent ou se modifient sans répit. Tout l'édifice conceptuel, rationnel qu'est le texte juridique serait donc miné en son fondement par la métaphore. Mais une telle approche n'est-elle pas elle aussi exposée à la critique ? En effet, si les concepts sont originairement métaphoriques, est-ce que le travail du juriste ne consiste justement pas à limiter le glissement incessant des dénominations et des concepts pour parvenir à une détermination rationnelle ? Le discours juridique s'efforce de travailler à partir de concepts précis dont le sens doit être délimité afin de bannir les ambigüités et les fausses évidences qu'elles véhiculent. La démarche du juriste consiste alors à rapporter l'inconnu au connu, ce qui revient finalement à faire usage d'une métaphore, à désigner l'inédit par une dénomination déjà connue et à déplacer ainsi le sens attaché à cette dénomination, à le transporter.

Mais le problème est que le discours rationnel finit par prendre des métaphores pour des essences. Dans La Métaphore vive, dans sa huitième étude intitulée « Métaphore et discours philosophique », Paul Ricoeur souligne ainsi que ce n'est pas parce que la métaphore est à l'origine du concept qu'elle continue d'agir en lui. Si l'empreinte métaphorique est effacée sur la monnaie conceptuelle, alors il n'y a plus de glissement possible et d'effet inaperçu. Si le texte juridique est riche en métaphores « lexicalisées », une fois que cette opération de lexicalisation a eu lieu, une signification propre est produite qui devient indépendante par rapport au premier sens du mot. La question du rapport entre métaphore et texte juridique se noue finalement autour de l'idée d'analogie. Mais si ce jeu sur le semblable est possible, c'est parce qu'il y a du même, parce que l'identité fonde la ressemblance. Si ressemblance il y a entre diverses réalités, n'est-ce pas parce qu'on peut les référer en dernière instance à l'Idée ? Si nous reprenons l'image de abeillage dans le droit seigneurial, si ces différentes abeilles ont une ressemblance entre elles, c'est bien parce qu'elles ont à voir avec l'Idée d'abeille, avec une unité et une identité qui rendent raison de cette ressemblance. Si le « même » fonde le « semblable », si le concept est logiquement premier par rapport à la métaphore, celle-ci est peut-être déterminante dans l'ordre de la découverte. C'est ce que Ricoeur appelle la «métaphore vive» ou dans notre optique une métaphore 'in vivo'. Ce qui résulte de l'irruption d'une métaphore, ce n'est pas un savoir, un donné, mais une demande en concept. Cet écart sémantique en quoi consiste la métaphore sollicite la pensée. L'interprétation est à l'intersection de ces deux mouvances en ce qu'elle recherche la clarté du concept tout en se nourrissant de la force et de la richesse de la métaphore 'in vivo' qui relance la réflexion quand le concept l'a figée. Si la métaphore a une place décisive dans le texte juridique, c'est bien parce que la pensée spéculative se sait différente et capable d'assumer le défi que lui lance la métaphore qui demande toujours à être travaillée par le juriste en même temps qu'elle lui est nécessaire. Si le discours juriste a besoin des possibilités offertes par la métaphore, la rigueur du concept permet seule d'en éviter les effets indésirables.

Néanmoins, si la métaphore apparaît comme irréductible, en fait, son danger le plus grand n'est pas sa présence, mais le fait qu'on oublie sa présence, qu'elle s'efface comme métaphore et agit alors sans que l'on en puisse maîtriser les effets. Elle tend à identifier des registres hétérogènes et à introduire une confusion périlleuse. Ainsi, Le mot « fortune » dans l'expression « fortune de mer » ne se réfère pas à la richesse. Le mot est pris dans son sens latin d'« aléa ». Il s'agit d'un événement fortuit qui peut être heureux ou malheureux. Ici il est pris dans le sens d'évènement dommageable qui se produit au cours d'une expédition maritime tel, un naufrage. C'est donc un cas de force majeure propre au droit maritime.

Si les dangers de la métaphore en langue de spécialité sont réels, ne peut-on pas envisager un usage maîtrisé de celle-ci qui la subordonnerait au concept tout en servant sa cause. En ce sens, on peut mettre en évidence des fonctions de la métaphore dans le texte juridique. Elle a d'abord un rôle didactique en ce qu'elle permet de présenter de façon imagée, en se référant au sensible, des relations abstraites. La métaphore terminologique est aussi un moyen de tisser des liens entre des domaines apparemment distincts et de révéler la systématicité et les analogies présentes dans la réalité elle-même. Elle est un outil particulièrement efficace pour retrouver des identités de structure au-delà de la multiplicité des domaines concernés. La métaphore terminologique, enfin, permet au juriste de suppléer aux insuffisances du langage, lequel ne peut offrir de nouvelles dénominations qui correspondraient à une façon inédite d'appréhender la réalité. La démarche phénoménologique est intéressante à cet égard. Elle fait l'épreuve des limites du langage naturel une fois l'épochè opérée.

Ainsi, souvent, ce que l'on essaie de dire « n'existe jamais réellement » c'est pourquoi nous percevons la métaphore terminologique comme le réceptacle et pour ainsi dire la nourrice de tout ce qui naît. Pourtant, notre incapacité d'abstraire totalement, la finitude, la rendent omniprésente. La déficience est alors du côté du sujet parlant et écrivant qui a besoin de l'appui des images pour penser.

L'apparente aporie du changement linguistique

Nous l'avons compris, la langue de spécialité, elle aussi, est primitivement et foncièrement métaphorique. Ainsi, l'objectif dans le raisonnement par analogie est alors d'exploiter un appariement adéquat effectué entre les éléments constitutifs des domaines source et cible, de manière à transférer un certain nombre de connaissances du premier vers le deuxième.

Conclusion et perspective

Nous avons présenté deux types d'approches effectués en parallèle : une approche computationnelle pour la détection automatique des métaphores et pour l'aide à leur analyse, une approche linguistique comparative utilisant les mêmes corpus que ceux de l'étude computationnelle. Nous projetons maintenant une confrontation des approches afin d'asseoir la validité des modèles proposés et dans le but de créer un outil qui sera mis à disposition de la communauté scientifique.

  Le corpus littéraire contient environ 200 romans du XIXème siècle français ou portugais et leur traduction dans l'autre langue. Il a été constitué grâce aux textes disponibles sur les sites des Projets Gallica de la Bibliothèque nationale (BNF, http://gallica.bnf.fr/) et Gutenberg (http://www.gutenberg.org) et par numérisation des traductions ou des textes en langue originale correspondants. Le travail de numérisation a été effectué grâce à la plate-forme technologique de la Maison de l'orient et de la méditerranée de l'Université Lyon 2 (Patrick Desfarges, responsable).

  Sélection d'expressions métaphoriques et d'expressions non métaphoriques appariées en fréquence d'emploi De chacune de ces 3 bases (base de syntagmes nominaux extraits du corpus littéraire, base de syntagmes verbaux extraits du même corpus et base de syntagmes nominaux extraits du corpus journalistique), deux ensembles ont été sélectionnés. Le premier ensemble contient des expressions métaphoriques sélectionnées et validés par I. Oliveira. Le second ensemble a été obtenu par extraction automatique à partir de chacune des 3 bases et de façon aléatoire d'un groupe d'expressions non métaphoriques. Nous avons ensuite vérifié que chacune des paires de groupes (les expressions métaphoriques et non métaphoriques de chacun des corpus) sont appariées en fréquence d'emploi. Mesure du recouvrement thématique Enfin, pour chaque expression sélectionnée et pour chaque mot plein de cette expression nous avons calculé les 50 premiers mots associés donnés par le modèle ACOM 1 . Le programme que nous avons réalisé calcule ensuite le nombre de mots associés communs à l'ensemble des mots pleins de l'expression donnée 2 . Ainsi pour l'expression non métaphorique un visage rose l'ensemble de mots associés communs au mot visage et au mot rose comportent les mots blanc, rouge, main, fille, visage, rose et la mesure vaut donc 6. Nous avons ensuite procédé à une analyse statistique (ttest pour échantillons indépendants) comparative des résultats entre expressions métaphoriques et expressions non métaphoriques appariées en fréquence d'emploi. Les résultats sont significatifs pour chacune des bases et montrent donc un recouvrement entre mots associés nettement plus faible pour les expressions métaphoriques que pour les expressions non métaphoriques. La figure 1 donne les résultats pour les expressions nominales extraites du corpus littéraire. Nous donnons ici les valeurs du test sur les autres séries : (1) Moyennes du nombre de mots associés communs pour les expressions verbales métaphoriques M=0,806763 et non métaphoriques M=2,464213 du corpus littéraire (t-test(299)=-5,75359, p < 0,001), (2) Moyenne du nombre de mots associés communs pour les expressions nominales métaphoriques M=4,507576 et non métaphoriques M=8,871795 du corpus journalistique (ttest(81)=8,871795, p < 0,001).

Figure 1 iciFigure 1 :

 11 Figure 1 ici Figure 1: Moyennes du nombre de mots associés communs pour les expressions nominales métaphoriques M=1,903704 et non métaphoriques M=7,884615 du corpus littéraire. t-test(211)=-10.08, p < 0,001 Cette mesure du recouvrement thématique semble donc être un critère très prometteur pour le repérage automatique des métaphores.

o 1

 1 foule : arrivée, mille, feu, milliers, bras, capitale o tiède : eau, robinet o foule+tiède : bain, feu, fête, images, immense Ici, le mot bain commun à la fois aux expressions courantes bain de foule et bain tiède est détecté par le modèle comme un élément porteur des images qui émergent de l'expression métaphorique foule tiède.  La bataille de l'argent : o bataille : champ, issue, perdu, cheval, livrer o argent : armes, marchés, permet, problème, utilisation o bataille+argent : perdu, perdre, lutte, gagner, gagné. Ici les mots perdre, gagner et lutte sont proposés par le modèle comme des éléments porteurs de la cohésion sémantique de l'expression la bataille de l'argent et représentent effectivement ce qu'argent et bataille partagent au delà de leurs spécificités propres. Nous voulons dans l'avenir réaliser une expérience en grandeur nature afin de confronter les résultats du modèle à ceux donnés par une analyse linguistique et tester à la fois sa validité et son utilité pour les linguistes. La suite de l'article présente une analyse linguistique faite en parallèle à ce travail de modélisation et qui est donc un élément essentiel de cette confrontation programmée. La terminologie littéraire : le perspectivisme De notre corpus littéraire se profile non seulement la reconnaissance de la spécificité des langues littéraires dans leur diversité, mais aussi la reconnaissance de la spécificité de l'usage littéraire du langage. Un langage littéraire qui exprime le génie individuel de son auteur. Dans ce sens, la vérité dépend largement de la perspective choisie du réel. Il n'y a pas de réel en soi, pas de vérité absolue, pas de sens unique et fondamental qui épuiseraient les significations de l'être. Le choix de la perspective va dépendre de la sélection minutieuse des valeurs sémantiques que le lecteur aura privilégié dans son interprétation, dans sa perspective. Aucune perspective n'est purement logique, c'est-à-dire neutre, objective. Le discours littéraire se distingue ainsi par la richesse de ses images et le jeu sur les dénominations. En littérature, cette figure de rhétorique consistant en une translation, un transfert de sens risque à tout moment de rendre inopérant l'effort d'univocité et de rigueur. D'abord, la métaphore permet de désigner les opérations jusque là inaperçues de la conscience : la conscience « vise », « saisit », « suspend ». Toutefois, lorsqu'il s'agit de désigner la subjectivité absolue qui n'est aucunement constituée, la métaphore reste le seul recours possible. Cette subjectivité absolue est «quelque chose qu'il faut désigner métaphoriquement comme « flux », quelque chose qui jaillit « maintenant », en un point d'actualité, un point source originaire... » Le discours littéraire fait ici l'épreuve de ses limites, non parce qu'il n'est pas assez riche, mais parce que la subjectivité absolue n'est aucunement constituée et échappe au langage qui ne peut désigner que le constitué. Le constitutif lui échappe par nature. La métaphore ne devient véritablement intéressante que lorsqu'elle ne s'adapte pas parfaitement à ce qu'elle entend désigner et oblige parfois l'écrivain à commenter et interpréter le dispositif qu'il a lui-même élaboré par le biais d'une métaphore filée. Il prend appui sur elle en même temps qu'il en marque les limites et parvient ainsi à nous amener à une compréhension de la réalité qu'il désigne.

4.

  figure

la copie d'une décision de justice ou d'un acte notarié

  Les concepts ne sont que d'anciennes métaphores usées, éteintes voire mortes. Stabilisés, ces libres regroupements deviennent des catégories qui peuvent prendre la forme de collocations, locutions, de cooccurrences, d'expressions idiomatiques, de catachrèses.Déf: L'expression s'utilise en droit pénal pour signifier que le public est exclu de la salle d'audience où se déroulent les débats. En 'évolution de la langue juridique dotée d'une personnalité à part entière résulte de l'ensemble des influences contingentes qui ont façonné son visage actuel. Comme tout organisme vivant, la langue est le produit d'une histoire avec ses hasards et ses ajustements. Dans le domaine juridique, une étude en diachronie s'impose pour pouvoir examiner avec plus de fiabilité le cheminement des termes, l'interaction entre les changements linguistiques et les milieux socioculturels. À la lumière de cette analyse, il apparaît clairement que l'étymologie est une forme d'étude diachronique. Pour comprendre le mécanisme du raisonnement analogique et son rapport avec la similitude, on doit s'intéresser impérativement à l'évolution du terme juridique où l'on recherche des situations potentiellement appariables avec la situation à analyser.

	TERME: HUIS CLOS	INTERSECTION/LIEN ANALOGIQUE
		"Huis" vient de l'ancien français signifiant
		"porte" d'où le vocable "huisserie" pour
		désigner l' entourage des portes et des
	procédure civile on emploie l'expression "EN	fenêtres.
	CHAMBRE DU CONSEIL" et on dira qu'en
	TERME: FRUITS CIVILS matière gracieuse les débats ont lieu en	INTERSECTION/LIEN ANALOGIQUE
	Définition : Les fruits civils qui sont obtenus chambre du Conseil.
	grâce à un contrat dont le capital est l'objet,	PRODUCTION, RENDEMENT.
	tel les loyers et autres revenus en argent
	procurés par une chose. 2. Par extension,
	qualification donnée à des biens qui
	constituent des produits au sens exact (dans
	la mesure où ils entament la substance de la
	chose qui les fournit) mais dont la production
	comme celle des fruits est périodique et
	régulière.
		INTERSECTION/LIEN ANALOGIQUE
		Le nom de "grosse" proviendrait de ce qu'à
		l'époque où les documents de justice étaient
	comportant la formule exécutoire.	rédigés avec une plume d'oie, les commis
	Elle est apposée par le Greffier de la juridiction qui a rendu la décision ou par le notaire qui a dressé l'acte contenant une	des Greffiers et des Notaires étaient payés au ROLE, de sorte que leur rémunération
	reconnaissance de dette.	étaient d'autant plus élevée que la copie était
		longue. Ils avaient donc tout intérêt à écrire
		en grosses lettres. On parlait alors d'"
		écritures grossoyées".

TERME: LA GROSSE DEF. La "grosse" est le nom donné à L

Si l'expression contient plus de deux mots pleins le programme calcule la moyenne deux à deux du nombre de mots associés communs aux mots pleins de l'expression considérée..


F Gimbert

D Amitrano

J Weiss

Crossover from quasi-static to dense flow regime in compressed frictional granular media

Keywords: 62.20.-x -Mechanical properties of solids 64.60.De -Statistical mechanics of model systems (Ising model, Potts model, field-theory models, Monte Carlo techniques, etc.) 91.60.Ba -Elasticity, fracture, and flow

Being ubiquitous in a large variety of geomaterials, granular assemblies play a crucial role in the mechanical stability of engineering and geophysical structures. For these applications, an accurate knowledge of the processes at the origin of shear localization, i.e. faulting, in frictional granular assemblies submitted to compressive loading is needed. Here we tackle this problem by performing discrete-element numerical simulations. A thorough analysis of the evolution of multiscale mechanical properties as approaching sample macroscopic instability is performed. Spatial correlations operating within the shear stress and strain fields are analyzed by means of a coarsegraining analysis. The divergence of correlation lengths is reported on both shear stress and strain fields as approaching the transition to sample instability. We thus show that the crossover from a quasi-static regime where the sample deforms infinitely slowly to a dense flow regime, where inertial forces play a significant role, can be interpreted as a critical phase transition. At this transition, no shear band of characteristic thickness can be defined.

Introduction. -The mechanical behavior of granular materials is of wide concern, from natural hazard to engineering applications. However, the evolution of properties towards the flowing instability is still partially understood. For assemblies of non-frictional, hard (nondeformable) and spherical particles submitted to a macroscopic shear stress, the concept of jamming [1] provides a powerful framework to understand the onset of granular flow. Force chains, i.e. heterogeneous distributions of contact forces on a scale much larger than the typical particle size, allow assemblies of non-frictional particles to exhibit jammed states [2]. These jammed states correspond to assemblies resisting small stresses without irreversible deformation, i.e. to assemblies that exhibit a macroscopic friction inherited from cooperative effects, whereas unjammed states are associated to irreversible flow under any applied shear [3,4]. The jamming transition for such spheres is controlled by the critical packing fraction φ [1]. (a) Present address.

Here we investigate the transition of granular assemblies toward macroscopic instability in a different situation. Elastically deformable and frictional disks are loaded under multi-axial compression. This compressive loading situation is relevant when studying geophysical instabilities such as the mechanics of granular gouges within fault zones or the dynamics of landslides. In the 2-dimensional configuration we consider, multi-axial loading consists in increasing the axial stress σ 1 while keeping the radial stress σ 3 , i.e. the confining pressure, constant (see fig. 1). Doing so allows sample shearing, as the deviatoric stress τ = σ 1 -σ 3 is continuously increased.

The differences between this configuration and classical ones used to study jamming [1] are twofold. First, in assemblies of frictional grains, the parameter that controls whether the grain assembly behaves as a jammed or an unjammed state is the fraction of grains that carry forces, named non-rattler grains, rather than the packing fraction [5]. Second, instead of keeping the sample volume constant during mechanical testing, the compressive loading conditions considered here involve a confining pressure that prevents the non-rattler-grain fraction from evolving freely. As a consequence, a percolating strong force network remains in the flowing phase, called the dense flow regime [6].

To investigate the behavior of geomaterials at the approach of macroscopic instability under these compressive multi-axial conditions, laboratory experiments have been previously conducted on either continuous rocks [7,8] and discrete materials such as sand [9] or synthetic analogous media [10]. In continuous materials, the bifurcation theory has been widely used to describe the macroscopic instability. This theory considers a transition from homogeneous to heterogeneous deformation materialized by the creation of a perennial macroscopic shear band of characteristic thickness, spanning the whole sample [11,12]. Within granular materials, either in experiments [9] or in simulations [13], deformation fields computed over large macroscopic strain windows in the vicinity of macroscopic instability also show perennial shear bands, with a possible characteristic size. On the other hand, by considering relatively small temporal scales under similar loading configurations, other studies report heterogeneous and long-range correlated kinematics of quasi-static granular flow [14,15]. At which temporal and spatial scales and at which stage of the loading does the granular assembly deform homogeneously? What are the relevant key features of the stress and strain fields associated to the onset of macroscopic instability? Here we investigate this problem using numerical simulations. We thoroughly quantify the heterogeneity of stress and strain fields at the approach of the transition from a quasi-static regime, where the sample resists to the applied stress by deforming infinitely slowly, to a dense flow regime, where inertia comes into play.

Simulation setting. -The numerical simulations we conduct use the molecular-dynamics discrete element method [START_REF] Radjai | Discrete Numerical Modeling of Granular Materials[END_REF]. Two-dimensional granular assemblies of a number N g of frictional circular grains are considered. To characterize sample size effects [17], 320 simulations with N g = 2500, 80 simulations with N g = 10000 and 20 simulations with N g = 45000 have been performed. The results presented here concern 10000 grains samples. The grains areas are uniformly distributed. The largest grain diameter D max is set such that D max =3 D min . Grains interact via linear elastic laws and Coulomb friction when they are in contact [START_REF] Cundall | [END_REF]. The normal contact force f n is proportional to contact interpenetration δ via f n = k n × δ,w h e r ek n is the normal contact stiffness coefficient. The tangential component f t of the contact force is proportional to the tangential elastic relative displacement, with a tangential stiffness coefficient k t = k n . Neither cohesion nor rolling resistance is considered at grain contacts. The Coulomb condition |f t |≤µ micro f n , where µ micro is the grain friction coefficient, requires an incremental evaluation of f t every time step, which leads to some amount of slip each time one of the equalities f t = ±µ micro f n is reached. A normal viscous component opposing the relative normal motion of any pair of grains in contact is also added to the elastic force f n to obtain a damping of the dynamics.

The simulation procedure is divided into two steps. First, the granular assembly is built by performing an isotropic compression of dilute frictionless grains sets. This method builds dense and highly coordinated initial packings, characterized by a packing density of φ i ≈ 0.85 and a backbone coordination number of z * i =4 ,w h e r e z * i is expressed as [19], N c being the total number of contacts and x 0 the fraction of rattlers grains. Second, multi-axial compression tests are conducted using µ micro =1 . W ep e r f o r mt h ee x t e rnal mechanical loading on the grain assembly by using periodic boundary conditions [START_REF] Radjai | Discrete Numerical Modeling of Granular Materials[END_REF]. A rectangular, periodic, simulation cell is defined by writing the linear operator h as h = h x 0 0 h y ,w h e r eh x and h y correspond to the size of the cell period in the radial and axial direction (see fig. 1). Then, stresses are applied solving for the dynamic equations of motion of h by ensuring that the internal stresses computed over the whole grain assembly following eq. (1) (see below) are equal to the prescribed external stresses σ 1 and σ 3 . We increase the axial stress σ 1 at constant rate by imposing a stress increment δσ tr 1 at each discretization time interval t r = mmin kn /25, where m min is the mass of the lightest grain. This stress control loading mode avoids stress relaxations and associated feedbacks that would be obtained under straincontrolled loading, i.e. adjusting δσ tr 1 order to axially deform at constant rate ǫ1 ,w h e r eǫ 1 is the axial deformation. In this stress-controlled case, stress and strain localization structures develop freely. We keep the confining pressure σ 3 as constant and sized by setting the contact stiffness κ = k n /σ 3 equal to 1000 [START_REF] Radjai | Discrete Numerical Modeling of Granular Materials[END_REF]. While chosen relatively low to ensure moderate computational times, this value of contact stiffness remains relevant for applications in geophysics and engineering mechanics. Within the Earth's crust, wave speeds [20] suggest κ-values of approximately 10 3 from hundreds of meters to several kilometers depth. Moreover, compression experiments performed on assemblies of glass beads of approximate Young's modulus E = 70 GPa submitted to 100 kPa of confining pressure lead to a κ-value of 700.

z * i =2 N c /(N g (1 -x 0 ))
Results. -

Macroscopic behaviour.

Figure 2 shows the macroscopic response of a granular sample loaded using δσ tr 1 =1.10 -6 σ 3 . The inertial number I quantifies the ratio between inertial forces and imposed forces and is defined as I =ǫ 1 m/σ 3 [START_REF] Radjai | Discrete Numerical Modeling of Granular Materials[END_REF], where m is the average grain mass. Initially, I is of the order of 10 -6 -10 -5 (see fig. 2(a)). As τ increases towards 2σ 3 , while undergoing fluctuations associated to plastic events, I remains lower than 10 -4 , which is often considered an upper bound for quasi-static conditions [START_REF] Radjai | Discrete Numerical Modeling of Granular Materials[END_REF]21]. Hence, in the range 0 ≤ τ<2σ 3 , the sample undergoes quasi-static deformation. At values of τ larger than τ c ≈ 2σ 3 , I increases by several orders of magnitudes, reaching values of about 10 -3 -10 -2 . This drastic change in I-values defines the macroscopic instability and the transition to a dense flow regime, where inertial forces become significant. , where 0 and 1 denote the initial and current configuration. Sample contraction is observed in the early stages of mechanical testing and essentially results from contacts elasticity, as it would not be observed in the limit of rigid grains [START_REF] Combe | i nDeformation Characteristics of Geomaterials, 3me Symposium sur le Comportement des Sols et des Roches Tendres[END_REF]. Sample dilation irreversibly occurs after the maximum contraction and the dense flow transition observed at τ c ≈ 2σ 3 occurs at packing fractions larger than the initial one.

The critical shear stress value τ c ≈ 2σ 3 materializes the softening of the granular assembly, as it corresponds to a macroscopic friction coefficient of µ macro ≈ 0.5t h a ti s significantly smaller than the friction µ micro = 1 considered at the grain scale. While macroscopic hardening can be associated to the non-zero macroscopic friction obtained for an assembly of frictionless circular particles submitted to a similar compressional loading [3], the assemblies of frictional particles herein considered show a reverse behaviour, characterized by macroscopic softening. In the present configuration, collective effects and structural rearrangements generate weak zones whose spatial extent that are much larger than the particle size grows towards the macroscopic instability. Inelastic deformation and flow concentrate within these weak zones, resulting in a macroscopic softening. The quantification of these effects towards the macroscopic instability through the analysis of stress and strain fields is the backbone of this manuscript.

Multi-scale analysis. We first characterize the stress field by quantifying the spatial extent of stress concentration regions using a coarse-graining analysis [START_REF] Marsan | [END_REF]24]. An average shear stress rate τ (see footnote 1 ) is computed at different stages of mechanical testing (cf. color dots in fig. 2) over a time window T and over a broad range of spatial scales L, ranging from the micro-scale corresponding to mesh size (see fig. 3 (left)), to the macro-scale corresponding to sample size. Subsystems of the granular assembly are selected by defining square boxes of size W (see fig. 3 (right)). Grains are considered to belong to a given subsystem if their mass centers are included in the box of size W . For a given box of size W , the actual area spanned by the mesh associated to the selected grains assembly slightly differs from the area of the box W 2 .T h i s difference is due to the irregularity of the contour of the mesh at the boundary of the domain. Thus, for a set of N box boxes of size W , we define the average scale L as

L = 1 N box N box k=1 L k ,
w h ereL k is the scale associated to the box number k, computed as the square root of the sum of mesh element surfaces. The smallest subsystem of average size L =1.1D max shows a maximum standard deviation of L k -values of 0.2D max for the Voronoi triangulation used to compute stresses and of 0.13D max for the Delaunay triangulation used to compute strains (see fig. 3). No overlap of scales L k occurs through successive values of W .F o r a given assembly of grains included in the box number k, the stress tensor is computed over all contacts carried by each grain using [25], i.e. writing 2)). We find ντ =1 .3, C =0 .5a n dδ =1 . B o ttom: τ vs. L at ∆ = 0.005 when increasing the timescale t from t = T to t =1 6 T. For graphical convenience, all computed values have been normalized by the ones computed at the micro-scale.

σ k ij = 1 L 2 k g k c k (r c k i -r g k i )f g k c k j , (1)
where f g k c k j is the j-th component of the contact force exerted on grain g k at contact c k , r c k i is the i-th component of the position vector of c k ,a n dr g k i is the i-th component of the position vector of the mass center of grain g k .T h e stress rate tensor is obtained by differentiating the stress tensor components calculated at two successive configurations. These successive configurations are selected in time and separated one from the other by a time interval T = N g × 100 × t r that corresponds to the travel time of elastic waves through the whole sample [17].

Figure 4 (top) shows the results of the multiscale analysis performed on the stress rate field from the early stages of biaxial testing up to τ c ≈ 2σ 3 . Curves are selected with respect to the control parameter ∆ defined as ∆ = τc-τ τc . At the early stages of macroscopic deformation, i.e. at large values of ∆, and at small spatial scales, τ decreases as L increases. Moreover, for L-values larger than a crossover scale l * τ , τ remains constant. We conclude that associated shear stress rate fields are heterogeneous for l ≪ l * τ and homogeneous for l ≫ l * τ . Associated fields provided on the right-hand side of fig. 4 illustrate the connection between l * τ and the size of spatial zones of stress concentration. For ∆ = 0.42, the stress rate field shows a roughly homogeneous pattern at large spatial scales, i.e. scales larger than l * τ ≈ 3D max in that case, while the region of stress localization observed on the top-left part of the snapshot exhibits a spatial extent of the order of l * τ ≈ 3D max that corresponds to the crossover scale of the associated curve. We interpret l * τ as the correlation length [24]. As macroscopic deformation proceeds, l * τ grows until reaching the size of the system at τ c .A tt h a t point, a power law scaling τ ∼L -ρτ is observed over the entire scale range, with ρ τ =0.38. We verified using various sample sizes that the cut-off remaining on the scaling at ∆ → 0 is a finite-size effect [17]. At τ = τ c , large and strongly localized structures characterize the shear stress rate field (see snapshot at ∆ = 0.005). These scalings suggest a progressive structuring of the stress rate field as approaching the transition to the dense flow regime. We report an associated divergence of the correlation length l * τ by performing a collapse analysis (inset of fig. 4 (top)) and formulating that l * = l * τ diverges as

l * ∼ L δ s ∆ ν L δ s + C , ( 2
)
where L s is the square root of the sample area and ν = ν τ =1 .3 ± 0.1 is the exponent of divergence. The parameters δ and C characterize the finite-size effect [17].

The divergence of l * τ ∼ ∆ -ντ at ∆ → 0 is confirmed from a similar analysis performed on other moments τ q of the shear rate [17], revealing the multi-fractality of the shear rate field at the critical point. These particular features of the shear stress rate field are only observed at the specific timescale t = T , that corresponds to the travel time of elastic waves. At increasing timescales, i.e. for t>T, the departure from the power law scaling observed in fig. 4 (bottom) denotes a progressive homogenization of the shear stress rate field, materialized by the decrease of l * τ as t increases (see also the associated snapshots). We interpret the multi-scale properties of the shear stress rate field observed at t = T to be associated to stress concentration events occurring simultaneously to the propagation of elastic information throughout the sample, which operates by means of elastic waves. The loss of scaling properties reported beyond this time, i.e. for t>T , implies that the multi-scale localization structures observed at t = T have little memory, limited to the travel time of an elastic wave. In other words, the loss of scaling properties at t>T is explained by the superposition of several stress concentration events operating at t = T , but uncorrelated in time one with the successive other.

We then investigate whether similar features characterize the shear strain field. A Delaunay triangulation is performed on the grain centers of the non-rattler-grains set (fig. 3). The partial derivatives are computed at the mesh scale as ǫ ij =1 /2(∂u i /∂x j + ∂u j /∂x i ), where (u 1 ,u 2 ) and (x 1 ,x 2 ) are, respectively, the incremental displacements and spatial coordinates of grain centers. By averaging the partial derivatives, an average value of the shear strain rate is computed at increasing spatial scales L using the coarse-graining analysis exposed previously. Using the constant timescale t = T to compute displacement rates, the correlation length associated to the shear Δ 0.65 0.42 0.28 0.18 0.12 0.07 0.03 0.005 Fig. 5: (Color online) Multi-scale analysis performed on the incremental shear strain field δγ. A selection of corresponding fields is shown on the right side: an arbitrary color scale (not shown) has been chosen for each snapshot. Top: δγ vs. L for decreasing values of ∆ towards the critical point. The deformation scale δǫ1 = δǫp =1 .10 -5 is used to compute values of δγ . The inset displays data collapse with respect to ∆ (eq. ( 2)). We find νγ =1 .3, C =0 .5a n dδ = 1. Bottom: δγ vs. L at ∆ = 0.005 when increasing δǫ1 from δǫ1 = δǫp to δǫ1 =3 2 δǫp. The inset displays data collapse with respect to ξ = δǫ 1 -δǫp δǫp (eq. ( 3)). We find αγ =0 .4. For graphical convenience, all computed values have been normalized by the ones computed at the micro-scale. strain rate field does not diverge as approaching the transition to the dense flow regime (not shown). Thus, the structure of the strain rate field does not form simultaneously in time with the stress rate field, which would be the case if one considered only the elastic component of the strain. Here, for N g = 10000, approximately 10 4 stress increments are prescribed during the propagation time of an elastic wave throughout the sample. Hence, a multitude of contacts, in our case about 5% of the whole contact network, are sliding although elastic information does not have time to travel across the entire sample. However, a progressive structuring of the shear strain field is observed when computing the scaling of the incremental shear strain δγ for constant macroscopic deformation windows δǫ 1 = δǫ p =1 .10 -5 . The divergence of the correlation length l * γ as approaching the transition to the dense flow regime is reported on the incremental shear strain field in fig. 5. This divergence is similar to the one observed on the shear stress rate field, as we find ν γ = ν τ =1 .3 (eq. ( 2)). For larger macroscopic deformation windows δǫ 1 >δǫ p , the deformation field no longer exhibits multi-scale properties and l * γ continuously decreases as δǫ 1 increases. This observation is in agreement with the study of [14], who reported a narrowing of the fluctuation velocity distributions at increasing timescales. As the material softens when τ i n c r e a s e s( s e efi g .2 ( a ) ) ,t h et i m eo f integration associated to the constant deformation window δǫ 1 = δǫ p decreases when approaching the dense flow transition. This time can be either smaller or larger than the travel time T of an elastic wave, depending on the imposed loading rate δσ tr 1 . All the presented scalings are indifferently obtained for loading rate values δσ tr 1 smaller than δσ tr 1 =1 .10 -6 . For larger values of δσ tr 1 , the divergence of the correlation length is no longer observed [17]. In that case, the mechanical behavior of the granular assembly is related to a dense flow regime from the initial stages of deformation, attesting that the scalings reported in this study are associated to the transition from quasistatic to dense flow regimes of deformation. Interestingly, δǫ 1 = δǫ p remains equal to 1.10 -5 whatever the loading rate considered, suggesting that a given, constant, amount of plastic activity has to operate in order to observe multiscale properties within the incremental shear strain field.

To better understand the role of plastic activity in setting the characteristic value δǫ p , we perform a four-point dynamic susceptibility χ 4 analysis on the inter-particles contact network [26,27]. From an "initial" contact configuration selected at an axial deformation ǫ init 1 , we compute the self-overlap order parameter Q ǫ init 1 (δǫ 1 )= 1 Nc Nc i=1 w i , where N c is the number of the non-sliding contacts of the initial configuration and w i is a step-function cutoff that equals 1 if no sliding event occurred at contact i over the deformation window ǫ init

1 → ǫ init 1 + δǫ 1 ,a n d0o t h e r - wise. The first two moments Q(δǫ 1 )= Q ǫ init 1 (δǫ 1 ) and χ 4 (δǫ 1 )=N c [ Q ǫ init 1 (δǫ 1 ) 2 -Q ǫ init 1 (δǫ 1 ) 2 ]o fQ ǫ init 1 (δǫ 1
) (calculated from sample-to-sample fluctuations) are computed for initial configurations selected at the previously used ∆-values, i.e. within the quasi-static region (τ<τ c ). This analysis evaluates the number and the associated spatial heterogeneity of sliding events that occur as the axial deformation increases. As no significant variations of Q(δǫ 1 )a n dχ 4 (δǫ 1 ) are observed with respect to ∆, fig. 6 shows the average functions Q(δǫ 1 ) and χ 4 (δǫ 1 ) computed over all the values of ∆. Major and minor force networks, defined by selecting contact forces larger and smaller than the average force, are considered separately. By construction, Q(δǫ 1 ) is initially equal to 1. As δǫ 1 increases, Q(δǫ 1 ) decreases but never reaches 0, as about 35% (respectively, 55%) of the contacts of the minor (respectively, major) network never slide. These non-sliding contacts associated to the remaining rigid bodies illustrate the localization of permanent deformation. At low values of δǫ 1 , χ 4 (δǫ 1 ) increases as χ 4 (δǫ 1 ) ∼δǫ β 1 with β =1.8. This increase results from the nucleation of spatially correlated events of contact sliding. As δǫ 1 exceeds the threshold value δǫ p =1 .10 -5 , χ 4 (δǫ 1 ) saturates. Thus, at δǫ p =1 .10 -5 , i.e. when multi-scale properties are observed on the shear deformation field (see fig. 5), all the spatially correlated contacts located close to the Coulomb criterion, i.e. susceptible to slide, have been destabilized. At this stage, only 4% (respectively, 3%) of contacts have slid in the minor (respectively major) network.

Conclusion. -We have shown that, in dense granular assemblies, stress and strain fields are both characterized by a growing correlation length that diverges as approaching the onset of macroscopic instability. This instability can therefore be identified as a critical point. A similar behavior has recently been reported in compressive failure of continuous materials [24,28]. Herein, this behavior is observed as long as the macroscopic instability is reached in the quasi-static regime. On the reverse, the critical behaviour is no longer observed for imposed loading rates large enough so that the inertial number is larger than 10 -4 before the instability [17]. We interpret these stress and strain specific structures as resulting from dynamic stress redistributions induced by the local dissipation of elastic energy, materialized by contact sliding. These features can only be observed for characteristic timescales for stresses, and characteristic macroscopic strain increments for strains. These characteristic timescales may be sensitive to varying initial packing properties, such as, for example, low coordinated and/or loose initial samples. This sensitivity has not been investigated in this study.

The last question that arises is whether a limit in decreasing correlation length l * γ is reached at ∆ → 0f o r values of δǫ 1 much larger than δǫ p . If so, this limit value observed on the shear strain field would characterize the thickness of a perennial macroscopic shear band potentially formed at the onset of instability. We hypothesize that, close to the critical point (∆ → 0), l * γ varies as

l * γ ∼ L δ s L δ s Π αγ + C , (3)
where Π = δǫ1-δǫp δǫp and α γ is the exponent of divergence with respect to Π. This hypothesis is tested from a collapse analysis (inset of fig. 5). We find α γ =0 .4, meaning that l * γ keeps decreasing as the considered deformation window size δǫ 1 increases. Thus, no intrinsic scale of saturation, potentially associated to a shear band thickness, is identified at the onset of macroscopic instability. * * *

1 Fig. 1 :

 11 Fig. 1: (Color online) Multi-axial loading configuration on a sample made of 225 grains. Empty circles correspond to sample replication.

Fig. 2 :

 2 Fig. 2: (Color online) Macroscopic response of a 10000 grains sample compressed at constant stress rate δσ tr 1 =1 .10 -6 σ3. The inertial number I and the sample surface variation ∆S/S0 are defined in the main text. Color dots correspond to color lines in figs. 4 (top) and 5 (top). Dashed lines indicate the limit between quasi-static and dense flow deformation regimes.

Figure 2 (

 2 b) shows the surface change of the granular assembly ∆S/S 0 = (hxhy)1-(hxhy)0 (hxhy)0

Fig. 3 :

 3 Fig. 3: (Color online) Left: Delaunay (top) and modified Voronoi (bottom) tesselations for a polydisperse granular material. Right: coarse-graining analysis on a 2500 grains sample.

Fig. 4 :

 4 Fig. 4: (Color online) Multi-scale analysis performed on the shear stress rate field τ . A selection of corresponding fields is shown on the right side: an arbitrary color scale (not shown) has been chosen for each snapshot. Top: τ vs. L for decreasing values of ∆. Configuration locations on the stress-strain curve are shown in fig. 2. The timescale T is used to compute values of τ . The inset displays data collapse with respect to ∆ (eq. (2)). We find ντ =1 .3, C =0 .5a n dδ =1 . B o ttom: τ vs. L at ∆ = 0.005 when increasing the timescale t from t = T to t =1 6 T. For graphical convenience, all computed values have been normalized by the ones computed at the micro-scale.

Fig. 6 :

 6 Fig. 6: (Color online) Susceptibility analysis performed in the quasi-static region on the sliding contacts belonging to the major (red curves) and minor (black curves) network. The vertical dashed line indicates the deformation value δǫp =1.10 -5 .

The invariant x,w h e r ex represents τ for the shear stress or γ for the shear strain, is computed as x = x Ix II ,w h e r ex I and x II are the principal components of the stress or strain tensor.

We thank Gaël Combe for providing the DEM code and for fruitful discussions. We thank Jean Braun for providing routines to compute Voronoi tesselations. All computations were performed at SCCI-CIMENT Grenoble.


You Zheng
email: you.zheng@utbm.fr

Oumaya Baala
email: oumaya.baala@utbm.fr

Alexandre Caminada
email: alexandre.caminada@utbm.fr

Alexandre Gondran
email: alexandre.gondran@enac.fr

A New Model For Indoor WLAN Positioning System

Keywords: combinatorial optimization, WLAN planning, positioning system

Introduction

In this paper we attempt to answer the following question : how to deploy a WLAN in order to guarantee the requested Quality of Service (QoS) while reducing the location error ? Such a problem includes two aspects : WLAN planning and positioning error reduction. To provide users a optimal wireless access to their local network, WLAN planning not only consists in selecting a location for each transmitter and setting the parameters of all sites, but also acts on allocating one of the available frequencies to each Access Point (AP) configuration [START_REF] Baala | The Impact of AP Placement in WLAN-Based Indoor Positioning System[END_REF]. And toward the indoor positioning system, once the Received Signal Strengths (RSSs) from all visible APs are measured and inputted, the location is estimated and outputted using the RSS distribution and machine learning technique [START_REF] Gondran | Wireless LAN planning : a didactical model to optimise the cost and effective payback[END_REF]. We propose a new approach where WLAN planning and positioning error reduction are modeled as an optimization problem and tackled together during WLAN planning process.

Modeling the Problem

AP model : Since different types of APs have different parameter values, we predefine a list of AP types provided for user choice. When one AP is selected, we define its set of parameters such as, azimuth, emitted power and frequency. Furthermore, we predefine a finite set of candidate sites where an AP may be assigned.

Received Signal Strength model : RSS is the principle parameter used in bit rate calculation and positioning estimation. We define three kinds of RSS thresholds. The interference threshold is the minimum RSS a client can receive as interference. The positioning threshold is the minimum RSS a client must received for positioning estimation. The reception threshold is the minimum RSS a client must received to allow it to communicate.

Traffic model :

The traffic model of the network defines how to represent the demand for network load it also provides a framework for expressing the QoS desired. In this model, we define service zone, which are represented by polygons covering parts of the building. Several characters of service zone are defined : number of users inside the service zone, bit rate in kbps desired by user for the service zone. In Our model the calculation area is decomposed into a grid. Each grid node is called Marking Position (MP). Thus these service zones are represented by the MPs inside these service zones which are called Test Points (TPs). In the same way as for the service zones, we defined forbidden zone, which must not be covered.

Positioning model : As the service zone defined in the traffic model, in the positioning model we define positioning zone where the positioning error reduction is considered. And these positioning zones are also represented by the MPs which are called Reference Points (RPs). For each RP, we find a set of MPs called Aliasing set in which all the MPs have the same RSS vector as RP.

Throughput model : For the network demand estimation, the QoS is evaluated by an indicator called Signal to Interference plus Noise Ratio (SINR) [START_REF] Baala | The Impact of AP Placement in WLAN-Based Indoor Positioning System[END_REF]. The definition is local for each user, given by : SINR = P Best RSS P Other RSS × γ(∆f ) + N Where : P Best RSS is the highest Received Signal Strength (RSS). P Other RSS are the RSSs of interference channels. γ(.) is the protection factor corresponding to the attenuation coefficient between channels. For the positioning error estimation, the accuracy of positioning is estimated by an indicator called Refined Specific Error Ratio (RSER) [START_REF] Gondran | Wireless LAN planning : a didactical model to optimise the cost and effective payback[END_REF]. The expression of RSER at a MP is as follows :

RSER = 1 (n(n -1)) i j =i dist(i, j)
Where : i, j : The sequence number of the MPs in the set Aliasing. n : The number of MPs in the set Aliasing. dist(i, j) : The Euclidean distance between the MPi and the MPj.

Optimizing the problem

Our optimization problem is to determine the decision variables (sites, the transmission powers, azimuths, frequencies) in order to minimize the financial cost under some constraints on QoS and maximal cost installation threshold. The financial cost we consider as objective includes the cost of the network deployment and the cost of penalty for unsatisfied demands [START_REF] Baala | The Impact of AP Placement in WLAN-Based Indoor Positioning System[END_REF]. Mathematically, the objective function is formulated by :

Site C Site network + T P β(p tp ) × ∆ tp + RP γ(p rp ) × E rp
Where, C Site network is the network installation cost of a site. β(p tp ) is the penalty coefficient according to the priority assigned to the TP. ∆ tp is the deviation between the require bit rate and the real bit rate. γ(p rp ) is the penalty coefficient according to the priority assigned to the RP. E rp is the magnitude of positioning error. In this function, the first term is the network installation cost, the second term is the cost of unsatisfied demands on QoS and the third term is the cost of unsatisfied demands on positioning. To solve this optimization problem, a Meta-Heuristics algorithm is being implemented.


Clément Peyronne
email: clement.peyronne@capgemini.com

Daniel Delahaye
email: daniel.delahaye@recherche.enac.fr

Laurent Lapasset
email: laurent.lapasset@capgemini.com

Marcel Mongeau
email: marcel.mongeau@recherche.enac.fr

Keywords: Conflict Resolution, Air Traffic Management, Optimization

Introduction

The most critical point of air-traffic control is to ensure safety separation between airplanes. To achieve this goal, a safety standard separation distance has been defined : 5 Nm (Nautical miles) horizontally and 1000 feet vertically. Air traffic controllers are responsible for ensuring the respect of these separation rules. Air traffic being constantly increasing, controllers in charge of an aviation sector must handle more and more flights. Nowadays, air-traffic management (ATM) has already used every available ressource to increase airspace capacity. From now to 2030, air-traffic is expected to be further multiplied by a factor two or three [START_REF] Alliot | L'ATM : 20 ans d'efforts et perspectives[END_REF]. Consequently, ATM will have to deal with this overload while ensuring at least equivalent standards of safety [START_REF] Alliot | L'ATM : 20 ans d'efforts et perspectives[END_REF].

Recently, two projects have been launched to find solutions for the future of air-traffic control : SESAR (Single European Sky ATM Research) in Europe, and NextGen in the United States. These projects aim at including more automation in ATM, introducing new technologies. Both SESAR and NextGen consider a new ATM concept where the objective is to keep each airplane as close to its original trajectory as possible. The SESAR European project aims at finding solutions to this problem by several actions on different time frames. One of the solutions SESAR encourages is to introduce more automation in air-traffic control, creating an automatic conflict resolution tool to provide a decision support to the controller. In this paper, we choose to begin with full automation as a first step, delaying the problem of providing conflict-free solutions to the controller to a later stage.

In the next section of this paper, we present our air-traffic conflict resolution problem. We describe our B-spline trajectory model in Section 3. Then, we introduce in Section 4 the new optimization formulation of our problem and the semi-infinite approach to tackle it. In Section 5, we detail the necessary reformulation to use derivative-free optimization method. After that, in Section 6, numerical results are presented and discussed and compared with previous results obtained using a stochastic optimization method (genetic algorithms denoted GA) presented in [START_REF] Peyronne | Optimizing B-splines using genetic algorithms applied to air-traffic conflict resolution[END_REF]. We draw conclusions and present our future work in Section 7.

The air-traffic conflict resolution problem

Our air-traffic conflict resolution problem is defined as follows. We are given the number of aircraft, noted N , involved in the conflict, and for each aircraft, the start and end points are given. We aim at calculating N optimal conflict-free trajectories. Optimality can be defined in several ways : for example, optimality in distance or in fuel consumption. Nevertheless, the particular choice of an objective function does not change essentially the way to tackle the problem. Remark that our previous stochastic methodology [START_REF] Peyronne | Optimizing B-splines using genetic algorithms applied to air-traffic conflict resolution[END_REF] attempted at minimizing the number of conflicts detected in a discretized space. All aircraft are assumed to start from their starting point at time t = 0. In this study, we work in two dimensions as we assume that each aircraft is constrained to follow its given altitude profile. In this study, we aim at solving conflicts considering a 20-minute time horizon and en-route aircraft flying at approximately 400 miles per hour.

Trajectory modeling

In this section, we introduce a B-spline trajectory model that relies on only one free control point per trajectory (more details in our previous paper [START_REF] Peyronne | Optimizing B-splines using genetic algorithms applied to air-traffic conflict resolution[END_REF]). The B-spline approximation is to find a curve approximating a set of points of R 2 called control points which completely defines the curve [START_REF] Guillod | Interpolations, courbes de Bézier et B-Splines[END_REF]. We choose to rely on B-spline for modeling trajectories because it is a very efficient tool for curve approximation in terms of both approximation quality and computational time. Moreover, B-splines feature interesting properties such as C 2 -continuity (crucial for modeling smooth aircraft trajectories).

Basically, B-splines are parameterized curves generalizing Bezier curve concept [START_REF] Peyronne | Optimizing B-splines using genetic algorithms applied to air-traffic conflict resolution[END_REF]. We define a parameter s ∈ [0, n] varying from 0 to n, n being the number of control points of the curve. In our model, we choose n = 3 control points : the start and end points (which are given) and a third control point which will be determined by our optimization variable vector. We can define the B-spline curve as the following 2D parameterized curve C :

C(s) = (σ x (s), σ y (s)), s ∈ [0, n]
Our optimization method will decide the location of the middle control point. In this case we choose to use only one control point subject to optimization for each trajectory (for each aircraft). We consider that more than one turn in this time window would be inconsistent with ATM requirements. The i-th aircraft trajectory is driven by a set of three control points called β i (u) where i = 1, . . . , N the aircraft index.

Semi-infinite programming formulation

In this section, we introduce a new optimization model for the optimal conflict resolution problem. Instead of going through this conflict-detection phase as in [START_REF] Peyronne | Optimizing B-splines using genetic algorithms applied to air-traffic conflict resolution[END_REF], the new idea here is to consider a continuous separation-distance constraint, making the approach much more direct. This constraint must be respected at all times, leading to an infinite number of constraints. Therefore, we propose modeling the problem with an objective function involving a finite number of optimization variable subject to an infinite number of constraints, giving rise thereby to a semi-infinite optimization problem as in [START_REF] Visweswariah | Noise considerations in circuit optimization[END_REF]. In this study, we focus on minimizing the total relative distance increase with respect to the direct route.

Recall that to model a trajectory we use B-splines involving one movable control point (the start and end points of the trajectory are not meant to move). Let u be the vector containing all these central control point locations (u i is the central control point of aircraft i). Let γ β i (s(t)) be the 2D-position of aircraft i considering its control point locations and time t, described by β i (u) (to clarify notation, from now on we will denote β i (u) as β i ).

Thus, our relative distance increase the objective function is :

f (u) = N i=1 T i (u i ) -L i L i
where T i (u i ) is the total travelled distance of the aircraft i following its modified trajectory and L i is the direct (straight line) route length from the start point to the end point.

Constraints must guarantee conflict-free situations. They have to ensure that all aircrafts are, at least, τ := 5 Nm away from each other (5 Nm is the standard separation norm in ATM). For two aircraft i and j, this constraints reads :

c ij (u; t) := γ β i (s(t)) -γ β j (s(t)) 2 ≥ τ , ∀t ∈ [0, t ij min ],
where ||•|| 2 is the l 2 norm, t ij min (u) = min(t i end (u), t j end (u)), and t i end (u) is the arrival time of aircraft i to its end point and is calculated from the B-spline trajectory defined through u. Therefore, we choose to formulate our problem as :

min u f (u) = N i=1 T i (u i ) -L i L i s.t. c ij (u; t) ≥ τ ,∀t ∈ [0, t ij min (u)]; i = 1, ..., N ; j = i + 1, N . ( 1
)

Reformulation for derivative-free optimization methods

Our problem clearly fits into the scope of semi-infinite programming (SIP), with time as the semi-infinite parameter. We first test one of the modern derivative free optimization (DFO) methods [START_REF] Conn | Introduction to derivative-free optimization[END_REF] in order to have a rough idea of the viability of our problem formulation. DFO do not attempt at approximating derivatives. More precisely, we choose a trust-region based DFO method using simple bounds and a minimum Frobenius norm model. For that purpose, we adapt our formulation to eliminate the constraints, since DFO methods are not designed to handle constraints (for now). We therefore include the constraints as a penalization in the objective function. In order to do that, we first reformulate the constraints to have a finite number of constraints In an analogous manner to [START_REF] Visweswariah | Noise considerations in circuit optimization[END_REF], we rewrite our constraint functions c ij (u; t) as :

c ij (u) = t ij min (u) 0 max{τ -c ij (u; t); 0} dt = 0.
We then include the new constraints in the objective function as a penalization, yielding to the new objective function f DFO as follows :

f DFO (u) = f (u) + ω N i=1 N j=i+1 c ij (u)
where ω is a user-defined positive weighting parameter. We shall see in the next section that the results we obtain are good.

Numerical results

In this section, we compare results obtained with the two optimization approaches (SIP and GA) on simple situations involving a very large number of aircraft. We perform all our tests on a 2.53GHz Intel(R) Core(TM)2 Duo on a Windows XP Operating System. We set ω = 0.01 for the weighting parameter that balances the trade-off between the distance and the conflict criteria. In this study, we consider simple, but difficult, academic conflict situations as test problems for our methodology. Each instance involves N aircraft uniformly distributed on a circle of radius 100 Nm in the first two situations, and 400 Nm for the last ones, where each aircraft is heading towards the center of the circle. We consider four test problems : N = 8, 16, 32 and 64.

In the following table, we display the objective-function value

f (u) = N i=1 T i (ui)-L i L i
obtained, the number of objective-function evaluations required (denoted fevals), and the calculation time (time) for both GA and the SIP approaches (the latter relying on DFO and denoted dfoSIP). Remark that, although the objective function is not exactly the same, we compare the same value : the total relative distance increase generated by the conflict resolution. In the above table, "F" denotes a failure (no conflict-free solution found) of the method within the allowed budget in terms of number of function evaluations (10000 for N ≤ 32 and 25000 for N = 64).

N
It appears clearly that the dfoSIP approach performs better than genetic algorithms in terms of solution quality and number of function evaluations. The calculation time however, is similar for both methods (inferior with GA on the 32 aircraft situation) . We also see that, as the size of the problem increases, the difference in quality between GA and SIP increases as well. One cause for this difference in performance, independent of the optimization method itself, is the way we implemented conflict detection in the genetic algorithms : the objective function used in GA (the number of conflicts computed via a discretization of space [START_REF] Peyronne | Optimizing B-splines using genetic algorithms applied to air-traffic conflict resolution[END_REF]) always overestimates the separation distance, making thereby the constraints harder to respect. This may be why GA does not find a solution for the 64-aircraft instance. To conclude, we remarked that GA and SIP yield similar solution on each of the smaller instances (8, 16 aircrafts) in terms of traffic organization with better length, CPU time and function-evaluation budget for SIP.

Perspectives and conclusion

We introduced a semi-infinite programming of the air-traffic conflict resolution problem that yields very encouraging results even though it does not seek for a global optimum. It appears that our simple dfoSIP approach already performs better than the GA methodology we developed in a previous study [START_REF] Peyronne | Optimizing B-splines using genetic algorithms applied to air-traffic conflict resolution[END_REF], where we also modeled trajectories via B-splines. B-splines allows us to reduce drastically the search space dimension.

We already have some work in progress on testing, less symmetric, test instances. We are currently implementing the objective and constraint analytic derivatives in order to be able to apply a classical local optimization method. We then intend to test both of our approaches on real air-traffic and especially on a complete day of traffic over France with a common way to compute the conflict detection (in other words, using exactly the same objective function for both GA and SIP).


Cyril Voyant

Gilles Notton

Christophe Paoli

Marie-Laure Nivet

Marc Muselli

Kahina Dahmani

Numerical weather prediction or stochastic modeling: an objective criterion of choice for the global radiation forecasting

Keywords: components; mutual information, stochastic, NWP, prediction

Numerous methods exist and were developed for global radiation forecasting. The two most popular types are the numerical weather predictions (NWP) and the predictions usingstochastic approaches.

We propose to compute a parameter noted  constructed in part from the mutual information which is a quantity that measures the mutual dependence of two variables. Both of these are calculated with the objective to establish the more relevant method between NWP and stochastic models concerning the current problem.

I. Introduction

Because of their random and intermittent trend, the integration rate of renewable energies is limited in the electrical grid. As a matter of fact, a limitation allows protecting the grid and to warrant a good supply Numerical weather prediction or stochastic modeling: an objective criterion of choice for the global radiation forecasting quality and to guaranty perfect production/consumption balance [START_REF] Mellit | Artificial intelligence techniques for sizing photovoltaic systems: A review[END_REF][START_REF] Paoli | Forecasting of preprocessed daily solar radiation time series using neural networks[END_REF][START_REF] Voyant | Optimization of an artificial neural network dedicated to the multivariate forecasting of daily global radiation[END_REF][START_REF] Voyant | Numerical weather prediction (NWP) and hybrid ARMA/ANN model to predict global radiation[END_REF]]. For increasing this insertion rate numerous solutions are studied and applied. Among these solutions, the main one consists in coupling the renewable energies with storage means (as hydrogen, battery, etc.). However, this coupling is not sufficient if the storage management is not well mastered. To do it, it is essential to anticipate the renewable energies production. Considering the grid manager's point of view (Köpken et al., 2004), needs in terms of prediction can be distinguished according to the considered horizon: following days, next day by hourly step, next hour and next few minutes.

With efficient prediction tools dedicated to grid managers, the PV part in the mix energy would be increased; actually in France, the intermittent energy contribution is limited to 30%.Several methods exist and have been developed for twenty years: numerical weather predictions (NWP) and predictions based on stochastic approaches [START_REF] Voyant | Numerical weather prediction (NWP) and hybrid ARMA/ANN model to predict global radiation[END_REF][START_REF] Sfetsos | Univariate and multivariate forecasting of hourly solar radiation with artificial intelligence techniques[END_REF]. If these two methods use mathematical approaches, the first one Voyant et al. models the atmosphere and oceans and allows predicting the weather from satellite images and primitive equations (nonlinear partial differential equations impossible to solve exactly through analytical methods) [START_REF] Paulescu | Weather modeling and forecasting of PV systems operation[END_REF]. The second ones consists in statistical models allowing generating alternative versions of the time series, representing what might happen over non-specific time-periods in the future [START_REF] Voyant | Optimization of an artificial neural network dedicated to the multivariate forecasting of daily global radiation[END_REF]. The choice between these two methodologies is not really scientifically established and the purpose of this paper is to propose an objective rule allowing guiding researchers or manager in their choices related to location and spatial and temporal resolution. This paper is organized as follows. The next section is dedicated to a brief description of the two types of global radiation prediction models, the NWP and the stochastic approaches. The section 3 introduces the methodology and the parameters proposed. Section 4 describes results and finally section 5 presents our conclusions and gives some perspectives to this work.

Numerical weather prediction or stochastic modeling: an objective criterion of choice for the global radiation forecasting

II. Prediction Models

The next part describes the two available types of models for the global radiation forecasting: the NWP and stochastic approaches.

a) NWP models

A numerical weather model is a computer program that simulates the atmospheric motion in space and time. A variety of weather phenomena can be analyzed and predicted by these types of NWP models. In this type of model, the atmosphere is represented by a 3D grid. The finer is the grid spacing the more elaborate is the simulation. The simulation by this model generates the future state of the atmosphere in each network points from its initial state [START_REF] Radnoti | The spectral limited area model ARPEGE/ALADIN[END_REF][START_REF] Voyant | Optimization of an artificial neural network dedicated to the multivariate forecasting of daily global radiation[END_REF][START_REF] Bouttier | The Météo-France NWP system: description, recent changes and plans[END_REF]Yessad, 2010]. Among all the NWP models, it may be mentioned weather

b) Stochastic models

The global radiation forecasting is the name given to the process used to predict the available amount of solar energy. Numerous predictive methods have been developed by experts around the world. Often the times series (TS) mathematical formalism is used [START_REF] Sfetsos | Univariate and multivariate forecasting of hourly solar radiation with artificial intelligence techniques[END_REF][START_REF] Mellit | Artificial intelligence techniques for sizing photovoltaic systems: A review[END_REF][START_REF] Paoli | Forecasting of preprocessed daily solar radiation time series using neural networks[END_REF][START_REF] Voyant | Optimization of an artificial neural network dedicated to the multivariate forecasting of daily global radiation[END_REF][START_REF] Voyant | Numerical weather prediction (NWP) and hybrid ARMA/ANN model to predict global radiation[END_REF]. It is described by sets of numbers that measures the status of some activity over time. It is a collection of time ordered observations x t , each one being recorded at a specific time t (period). A TS model ( t ) assumes that past patterns will occur in the future. TS prediction or TS forecasting takes an existing series of data x t-k , .. , x t-2 , x t-1 and forecasts the x t data values. The goal is to [START_REF] Mellit | Artificial intelligence techniques for sizing photovoltaic systems: A review[END_REF][START_REF] Paoli | Forecasting of preprocessed daily solar radiation time series using neural networks[END_REF]. Note that all the stochastic models are not equivalent, depending on the problem to solve (time horizon, spatial characteristics, locations, etc.) the best predictor to use can change. The absolute prioritization of these models is not possible.

x ˆx ˆ

c) Methodology

Our objective is to determine the rule optimizing the choice between NWP and stochastic models considering a given spatial and temporal resolution and a location. To reach this objective we propose to compute a parameter noted . It is constructed from the ratio between two subparameters, the first is related to the distances for which the global radiation series are independent () and the second one from a time lag for which there is no mutual dependence ().

    =      



We will see in the following how to compute these two sub-parameters, but at first let us consider We can highlight that the dimension of this parameter , is unconventional: pixel.time_lag -1 . The interpretation of this parameter should indicate if, for a given spatial and temporal

Numerical weather prediction or stochastic modeling: an objective criterion of choice for the global radiation forecasting resolution and for a given location, the TS formalism and the NWP are relevant. The most important parameter for stochastic models seems to be more it is important and more the information contained in the past series (intrinsic behavior) will be used to model the future. Concerning the global model NWP, it is mainly the spatial link between the points of the mesh which will guarantee the merits of the approach. In fact, the time dependence and the spatial dependence are linked, but the temporal aspect is certainly to a lesser extent. The kinetic of the cloud cover (related to the primitive equation and nonlinear partial differential equations) must be observed pixel by pixel in order to develop relevant model. Indeed, if the distance between mesh grid points is too high, local phenomena (cumulus cloud has a scale of less than 1 km) are not taken into account. But as all the pixels will be with the same average cloud cover, NWP will be very relevant in this case. For the high resolutions, the clouds local displacements appear: high heterogeneities between pixels are generated. Moreover, according to the dynamic and chaotic appearance of the fluid dynamics equations involved in weather forecasting, their modeling is very difficult and even impossible. Referring to the Equation (1), this spatial link will be illustrated by  Thus if stochastic models will be preferable and if >>1 NWP models will be more appropriate. We must note that considering the location, the time lag and the spatial resolution,  can fluctuate significantly. Mutual information is more general and measures the reduction of uncertainty in Y after observing X. So MI [START_REF] Kuijper | Mutual information aspects of scale space images[END_REF]

can measure

Numerical weather prediction or stochastic modeling: an objective criterion of choice for the global radiation forecasting non-monotonic and other more complicated relationships. It can be expressed as a combination of marginal and conditional entropies (respectively H(X) and

) Y H(X
) [START_REF] Lauret | Bayesian statistical analysis applied to solar radiation modelling[END_REF].

 ) Y X)-H(X MI(X,Y)=H(   



For the value of  we suggest computing and analyzing the mutual information between the global radiation and the distance between considered points and for the value of We will focus on the mutual information between the global radiation and the time lag for the considered location. Entropy corresponds to a measure of unpredictability or information content and can be written by the following expression (entropy of a discrete random variable X with possible values {x 1 ,..x n }):

    x x p x p X H )) ( log( ) ( ) (     Voyant et al.
One may also define the conditional entropy of two events X and Y (Equation ( 4)). This quantity should be understood as the amount of randomness of the random variable X given that you know the value of Y.

    x y y x p x p y x p Y X H )) , ( / ) ( log( ) , ( ) (  

 ,Y)/H(X) X,X) =MI(X I(X,Y)/MI( nMI(X,Y)=M

   

Numerical weather prediction or stochastic modeling: an objective criterion of choice for the global radiation forecasting

III. Results

The mathematical formalism is being enounced; the result concerning the spatial and temporal dependence of the global radiation in Corsica Island is exposed in this part. To illustrate the previous methodologies we decided to use in this study the HelioClim-3 database (HC-3). HelioClim is a family of databases containing solar irradiance and irradiation values at ground level. We have irradiance images (at an hourly step ; Wh/m²) over 8 years (2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012) for overall Corsica (1150 points of measurement, around 42°1'N and 9°E, with a pixel area lower than 6.5km², that is to say a grid spacing distance of 2.5 km) like seen in Figure 1 and Figure 2. We see that for a distance lower than 5 km, the nMI value is high, but it decreases after 5 km and it remains stable (~0.85%). The space parameter () is determined by the intersection between the limit of nMi and the tangent of the fitted curve at 0 (see Fig. 2). The chosen fitted curve is an exponential decay. In fact, the limit is close to the pixel size:

~ 2.5km. But for other countries (and so climates), for other time steps

 Voyant et al.
and for other mesh grid resolutions, nothing suggests that this 1 pixel limit will be maintained. Beyond this threshold, we can consider that the global radiations received on the other pixels are independent.

b) Temporal dependence

In Figure 4, we see for two points located in Ajaccio and in Corte, the MI versus the time lag (auto-mutual information). The first location is a seaside site and the second one is a mountainous site. In this example, we see that the first minimum is obtained for the 9 th and the10 th time lag, corresponding to 9 and 10 hours (as a complementary, for Corte the first minimum is 7 hours). Consequently using a stochastic model, for these localizations is inconsistent to predict the global irradiation at a moment t from data collected 9 and 10 hours before (in case of 1 hour time step and measured time series, see [START_REF] Voyant | Numerical weather prediction (NWP) and hybrid ARMA/ANN model to predict global radiation[END_REF]). If we compute the auto-mutual information for all the points of the meshing (1150 pixels), we see that the median value is 7 (min=5, max=10, mean value=7.63 and standard deviation=1.08). For the overall territory,  is around 7 time steps i.e. 7 hours.

c) Interpretation

In this study (overall territory) and with the methodology used, is time lags and is 1 pixel and consequently

1 7 / 1    .
In the Ajaccio and Corte cases, the  parameter is also less than 1 (respectively 0.1 and 0.11 pixel.time_lag -1 ). As said in the beginning of this paper (section 3), we can estimate that, in considering a hourly time-step and for a spatial resolution of 2.5 km, the use of a stochastic model is relevant. In the literature, this aspect is confirmed if the stochastic model gives a nRMSE close to 15% [START_REF] Sfetsos | Univariate and multivariate forecasting of hourly solar radiation with artificial intelligence techniques[END_REF][START_REF] Voyant | Optimization of an artificial neural network dedicated to the multivariate forecasting of daily global radiation[END_REF][START_REF] Voyant | Numerical weather prediction (NWP) and hybrid ARMA/ANN model to predict global radiation[END_REF][START_REF] Voyant | Numerical weather prediction (NWP) and hybrid ARMA/ANN model to predict global radiation[END_REF], NWP give less pertinent results (nRMSE> 20% [START_REF] Paulescu | Weather modeling and forecasting of PV systems operation[END_REF]). If the spatial resolution increases,  will also increase and  will be able to be close to one or more. If now we consider a prediction of the global radiation 24 hours ahead by hourly step for example, there are always two methodologies: the NWP or the stochastic approaches. The two methods are able to propose results more or less relevant, but to make a choice between the two forecaster types without test them (it is long and laborious), it is also possible to use the ad-hoc index defined in this paper (). We know that corresponds to time lags so the stochastic model will not be certainly relevant for the 24 hours horizons, however it no sure that the NWP model will be more efficient. In this case is yet Numerical weather prediction or stochastic modeling: an objective criterion of choice for the global radiation forecasting pixel. Note that in the 24 hours ahead case, the prediction during the night of the global radiation does not make sense, the night values can be deleted and the 24 hours prediction can be replaced be a 12 hours predictions and the global series by a new series where only the hours between 8AM (included) and 7PM are considered. Note that this approximation induces that during summer days, the hours of sunrise and sunset are not considered. According to the Figure 4, the maximum lag to considered is yet equal to 7 hours as for the 1 hour horizon (9 and 10 hours for Ajaccio and Corte points), the data are the same thus, as the horizon of prediction is equivalent to 12 hours, the time lag is equivalent to 0.6 time lags (~7/12). In this configuration, the index becomes higher than 1 (pixel.time_lag -1 for the overall territory, 1.2 pixel.time_lag -1 for Ajaccio and 1.3 pixel.time_lag -1 for Corte). The proposed interpretation induced that the stochastic models become not relevant, and the use of NWP is recommended. This  interpretation is verified in [Voyant, 2013] and in [START_REF] Paulescu | Weather modeling and forecasting of PV systems operation[END_REF]. In the Voyant et al. first reference the stochastic models (MLP and ARMA) give an nRMSE close to 30% in the h+24 case and the second anRMSE between 10 and 40% (average 25%) for the NWP model and this horizon. This methodology is relevant for other horizon, other mesh grids and other locations.

IV. Conclusion

In this paper, we proposed a methodology to justify the use of a NWP or a stochastic model according to three considered parameters: spatial resolution, temporal step and location. The mutual information and the three proposed parameters are the mathematical tools used as choice criterion between forecasting methodologies. For our case study (Corsica Island), we see that for a temporal resolution of 1 hour and a spatial resolution of 2.5km, a stochastic model is the best choice, but in the 24 hours head forecasting (prediction by hourly step for the next day) the Numerical weather prediction or stochastic modeling: an objective criterion of choice for the global radiation forecasting NWP forecasting based are the most relevant. In view to generalize this described methodology, we must validate it in other locations and for various spatial and temporal resolutions. Moreover a regional scheme will be performed, in order to separate the different microclimates and to develop model able to consider the very high resolutions of the nonhydrostatic models (as Meso-NH for example and the sub-kilometric meshgrid). In this case, a new formalism should be developed to take into account the fact that the spatial MI is non-continuous at the origin (nugget effect) and because the operated normalization (nMI) generates an offset of the curve. Another important improvement should be operated to clarify the threshold of , the limit fixed to 1 is not clearly proven: what happens for a value of  between 0.9 and 1.1? Is it really possible to conclude? It will be the objectives of a future paper.

VI. Conflits of interest

None

VII. Acknowledgement

Thanks to Armines and to the SoDa community for making the HelioClim data available (http://www.helioclim.org/index.html).

Numerical weather prediction or stochastic modeling: an objective criterion of choice for the global radiation forecasting

VIII. Figures

  research and forecasting model (WRF), AROME which concerns the mesoscale, or the model of the European centre for medium-range weather forecasts (ECMWF). The prediction errors of this model depend on the considered locations and fluctuate between 20 and 40% (nRMSE= normalized Root Mean Square Error, [6]).

  The question remains of how to find the values of and We propose to use the mutual information tool (MI(X,Y) in (2)) which is a quantity measuring the mutual dependence of two variables X and Y. In fact, this formalism replaces and generalizes the cross or auto-correlation and classical variogram concepts which allow to measures only the linear relationship between two variables X and Y (Pearson correlation and variogram) or two elements of time series X and L i X (L and ithe lag operator and associated order ; partial or normal autocorrelation factor).

   The definition of the joint probability distribution function (p(x,y)) and marginal probabilities (p(x) and p(y)), allows to define a new form of the mutual information as described in the equation (5). of the mutual information is often the bit (if the binary logarithm is used), thus it is possible to normalize this quantity by his maximal value to obtain a percentage. The new parameter is called the normalized mutual information (nMI, Equation (6)).

Figure 1 .

 1 Figure 1.Points of measurements meshgrid, circle locating Ajaccio (41°55'N and 8°44'E, elev. 0-787 m), square locating Corte (42°18'N and 9°09'E, elev. 300-2626 m) and triangle Bastia (42°42'N ; 9°27'E ; elev. 0 m)The global radiation time series used are not related to measurements

Figure 1 .

 1 Figure 1. Irradiance map computed during spring 2012 in Corsica (HC-3).The unit of the color map is Wh/m²

Figure 2 .

 2 Figure 2. Normalized mutual information of the global radiation versus the distance between considered points

Figure 3 .

 3 Figure 3. Mutual information of the hourly global irradiation versus the time lag for a 2 given points of the grid: Ajaccio (top) and Corte (bottom)

  Voyant C., Randimbivololona P., Nivet ML., Paoli C and Muselli M. 24 hours ahead global irradiation forecasting using Multi-Layer Perceptron.Meteorological applications, 2013 in pressYessad K., "Basics about ARPEGE/IFS, ALADIN and AROME in the cycle 37 of ARPEGE/IFS", Météo-France/CNRM/GMAP/ALGO, Research report, 2010

Figure 4 .

 4 Figure 4.Points of measurements meshgrid, circle locating Ajaccio (41°55'N and 8°44'E, elev. 0-787 m), square locating Corte (42°18'N and 9°09'E, elev. 300-2626 m) and triangle Bastia (42°42'N ; 9°27'E ; elev. 0 m)

Figure 5 .

 5 Figure 4.Points of measurements meshgrid, circle locating Ajaccio (41°55'N and 8°44'E, elev. 0-787 m), square locating Corte (42°18'N and 9°09'E, elev. 300-2626 m) and triangle Bastia (42°42'N ; 9°27'E ; elev. 0 m)

Figure 6 .Figure 7 .

 67 Figure 6. Normalized mutual information of the global radiation versus the distance between considered points